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ABSTRACT 


High energy laser weapons have been evolving progressively in recent years. These 
weapons deliver high-intensity beams to a target and can instantly destroy or burn it. 
They may cause potential threats to Navy ships, computer networks, guided missiles, and 
satellites in orbit. In order to reduce our military’s vulnerability to high energy laser 
weapons, one possible countermeasure is to rotate or rock the object itself when it is hit 


by the laser beam. 


The main purpose of this thesis is to investigate the relationship between the 
speed of a rotating/dithering laser beam and the maximum temperature rise induced by 
the laser beam on a finite solid. We have investigated extensively the numerical 
solutions for the transient temperature rise in both one-dimensional (1-D) and two- 
dimensional (2-D) finite solids due to rotating/dithering laser beams. Our mathematical 
approaches include the eigenfunction expansion method, the Crank-Nicolson method, the 
Fast Fourier Transform method, and COMSOL for 1-D and 2-D cases. We have 
employed COMSOL to solve the 3-D nonhomogeneous heat equation. 


This thesis provides the first study that we know of on the effect of 
rotating/dithering laser beams on a finite target. Our results are consistent with previous 
analytical studies on semi-infinite regions. The quantitative relationship between 
maximum temperature rise and laser beam rotating speed, which is presented in this 
thesis, can be used as a general guide for adjusting the speed of rotation of the target in 
order to prevent the maximum temperature rise from reaching the melting point of the 


target. 
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I. INTRODUCTION 


Laser-induced heating has been widely studied [1, 2, 3, 4]. Most of these studies 
are restricted to a moving/scanning laser beam. Recently, the effect of a 
rotating/dithering laser beam on the temperature rise of a semi-infinite domain has been 
studied [5]. However, in many realistic problems the geometry of an object is finite. So, 
in this thesis, we extend the previous studies to a finite domain and pinpoint the effect of 


a rotating/dithering laser beam on the temperature rise of a finite object. 


The laser beam can be considered as a moving disc heat source with a Gaussian 
distribution of heat intensity. The analytical solutions can be used to determine the 
temperature rise distribution in and around the laser beam heat source on the work 
surface as well as with respect to depth at all domains close to the heat source. The 
analytical and numerical solutions from mathematical approaches provide a better 
appreciation of the physical relationships between the relevant laser parameters and the 


thermal properties of the work piece. 


The main purpose of this thesis is to investigate the relationship between the 
rotating/dithering laser beam speed and the maximum temperature rise induced by the 
laser beam on a finite solid. In Chapters II and III, we present the analytical and 
numerical methods for obtaining a transient temperature distribution with various 
methods, which include the eigenfunction expansion method, the Crank-Nicolson method 
and the Fast Fourier Transform method (FFT). In Chapters III and IV, numerical 
solutions for one-dimensional (1-D) and two-dimensional (2-D) nonhomogeneous heat 
equations from MATALB and COMSOL are given. We compare different methods and 
investigate the relative errors of numerical solutions obtained in MATLAB and 
COMSOL. In Chapter V, we have employed COMSOL to solve the three-dimensional 
(3-D) problem. 


The objective of this thesis is to provide insights on developing objects of 
adequate rotating speed against direct energy weapons (DEWs). The simulation results, 


especially the quantitative relationship between maximum temperature rise and 


| 


rotating/dithering speed, can be used as a blueprint to adjust the speed of rotation of the 
target against DEWs to prevent the temperature rise from reaching the melting point of 


the target. 


I. ANALYTICAL SOLUTION FOR A TRANSIENT 
TEMPERATURE DISTRIBUTION DUE TO A ROTATING OR 
DITHERING LASER BEAM IN A FINITE SOLID: EIGEN- 
FUNCTION EXPANSION METHOD 


Consider a laser beam source used to heat a 1-D finite rod. A 1-D heat equation 


can be written as [1] 


Ou Ou a, 
— = g(t O0<x<L Eq. 2.1 
Cot ° Ox K, q02)) : aii 


where u(x,t) 1s the temperature rise with respect to the ambient temperature, g(x,t) 1s 
m 
the heat source created by laser beam, @, 1s the thermal diffusivity with units ae and 


S 
Sdn ._ | W 
K, 18 the thermal conductivity with units) | ; 
m 





Two assumptions are made in this analytical approach. First, heat losses by 
radiation are negligible as compared to the intensity of the incident laser beam. Second, 


thermal properties, such as thermal conductivity K, and thermal diffusivity @, , are 


considered constant and evaluated at an average temperature. The second assumption 


makes this heat equation a linear problem. 


For this 1-D rod, we impose insulating boundary conditions (BCs) 


(0,1) =(L,,1)=0 , which assumes that no energy escapes into the air at the 
x x 


air/material interface. This is a good approximation for most materials under 
consideration, as heat flow by conduction through the material far exceeds the heat loss 
by radiation or convection at the air/material interface. 

The initial condition (IC) is u(x,0)=0, which reflects the fact that there is no 


temperature rise with respect to the ambient temperature before the laser hits the rod. 


The heat source g(X,t) is considered as a laser beam of temporal continuous wave 
(CW) and spatially modeled by a Gaussian distribution, which corresponds to the 
theoretical JEM, mode of the laser. The term TEM comes from the acronym 


“transverse electromagnetic mode.” The subscript of TEM specifies the number of nodes 


generated by a slight misalignment of the mirrors located in the laser cavity [1, 6]. Figure 
1 shows some transverse modes and the simplest modes, TEM is used in this thesis; 


the heat source generated by the laser beam is a Gaussian distribution. 





TT Node 
TEM 5, 





Node 
Node 
.s$ Node 
TEM 
10 TEM, , 


Figure 1. | Different transversal modes in a laser spot [From 1 | 


For a dithering laser beam in JEM, mode, the heat source g(x,f) can be written 


as 


0-4, (0) oe 
q(x,t) = I,e : aS eae i= 


Here, x.(t)is the position of the dithering Gaussian beam and x, is the initial 
position of the laser beam. In the above formula, the center point of the rod is picked. 7 
is the effective radius of the laser beam. /, is the intensity of the laser beam at the center 


of the heat spot after it 1s absorbed by the material. k is a constant used for the Gaussian 


model [2]. Figure 2 depicts a dithering laser beam shining on a 1-D rod. 


k 2 
-—x- x, (4) 
Dithering Laser Beam g(x,t)=J,e ” 











L L L 
0 =O x — ag L 
Z 2 Z x 
T 
72 @t=0 Qe". 


Figure 2. Dithering laser beam on a 1-D rod 


We will apply the eignefunction expansion method to solve this 1-D heat equation. 
Generally speaking, the method of eigenfunction expansions consists of building up the 
solution of a boundary value problem as a sum of eigenfunctions of a Helmholtz problem. 
To begin with, consider a simple 1-D mathematical model in the Cartesian coordinate for 
solving the linear and homogeneous heat equation with initial conditions (ICs) and 


boundary Conditions (BCs) described above where no heat source g(x,t) is involved, the 


1-D heat equation becomes: 


eu =a V~u. 


Pa T 


To make the derivation process simpler, set Ap = 1. The method of separation of 
variables (SOV) is applied in solving this homogeneous heat equation [1]. We will move 
on to a nonhomogeneous solution with heat source g(x,t) by applying the method of 
eigenfunction expansion. 

It is intuitional that the temperature u is a function of position and time, 
therefore u = u(x,t). In the method of SOV, we attempt to determine solutions in the 


product form, 
u(x,t) = X (x)T(t) (Eq. 2.2) 


where X (xX) is a function of space x alone and T(t) is a function of time ¢ alone. This 
SOV reduces a partial differential equation (PDE) to several ordinary differential 


equations (ODEs). 


For a heat equation without heat source: 


SED = V~ u(x,t) (Eq. 2.3) 


We substitute the assumed product form of Eq. 2.2 into this heat equation and get 
X(x)P(t) = X"(x)P(t) 


We can actually “separate variables” by dividing both sides of this equation by 
X (x)T(t) to obtain 
Ss 


ae Ea. 2.4 
ae (Eq. 2.4) 


Now the variables have been “separated” in the sense that the left-hand side of the 
equation is only a function of timer and the right-hand side function is a function of 
space x. Since the variables t and x are independent of each other, the only way to get 


equality is to have the functions on both sides of (Eq. 2.4) constant and equal. Thus, 


TS Pn 
— —-—_ - -j 
Tr x 


where J is an arbitrary constant called the separation constant. Later we will explain 
why this separation constant is a negative value with A>0. Since we treat the laser beam 
as a heat source, and there is no heat source involved yet in solving the homogeneous 
equation, it is intuitive that the temperature decreases as a function of time, so 


PO __ 
a A<0 (Eq. 2.5) 


(Eq. 2.5) is a first-order linear homogeneous differential equation with constant 


coefficient. This ODE can be solved directly by seeking exponential solutions, 7’ = elt 


By substitution, the characteristic polynomial is r=—/. Therefore, the general solution 


of (Eq. 2.5) is 
T(t)= ce Ht (Eq. 2.6) 


Where Cc is an arbitrary constant. The time-dependent solution is a simple 
exponential. Recall that 7 is the separation constant, which for the moment is arbitrary; 
however, eventually we will find out that only certain values of A are allowable. 
IfA >0, the solution exponentially decays as ¢ increases because of the negative sign; 
ifA <0, the solution exponentially increases; if =0, the solution remains constant in 
time. Therefore, 7 > gives a reasonable solution since temperature decreases as time 


goes by. 
Our next move is to consider the right-hand side of (Eq. 2.4): 


xX " 
<-=-/ (A>0) (Eq. 2.7) 


This is a second-order, constant coefficient homogeneous ODE, substituting 


X (x) =e’ into (Eq. 2.7) yields the characteristic equation po LRAO, 


+iVax I 


Since A > 0, exponential solutions have imaginary exponents, X (x) =e n 


this case, the solution oscillates. Recalling Euler’s formula, el@ _ cos(@)+isin(@), so 
the choices COS (VAx] and sin(/Ax) are made to get the desired real number solutions. 
The general solution for (Eq. 2.7) 1s: 

X (x)= Acos(VAx)+ Bsin(VAx) (Eq. 2.8) 


This is an arbitrary linear combination of two independent solutions where A and 


B are just arbitrary constants. 


Both ends of the 1-D rod are insulated; the BCs on x direction states: 


Ou Ou 
ee =0 and — = 0 
OX |x=0 ox x=L,. 
This is equivalent to 
X\0)=0 and X\(L,)=0 (BCs) 


The derivative of (Eq. 2.8) is 

X (x) =—AVA sin(VAx)+ BVA cos(VAx) (Eq. 2.9) 

We can plug the first BC X'(0)=0 into (Eq. 2.9) to obtain the solution 
X (0) =0 which implies that B=0. So X(x) is a function of cosine only since the sine 
term vanishes. Another BC implies that Asin(VAL.)=0. To avoid the trivial 


solution X(0)=0, we take A=1, which forces sin(VAL.)=0. Since the sin 


function vanishes at the integer multiples of 7 , we conclude that 


Ja = JA, ==, n=0,41,42,....... 


Xx 


and so 


xX=X.(x%)= cos(—), 7 | ee eee 


Note that for negative values of n, we obtain the same solution; hence, solutions 


corresponding to negative n's may be discarded without loss. These X, (x) are from the 


eigenfunction of the 1-D heat equation. We now go back to (Eq. 2.6) and substitute 


2 
nit 
A =| — | to get 
a2 : 


_| nt 
L 


Xx 


t 








T(t) =c,e 


Now we return to (Eq. 2.1). We solve this problem using the method of 
eigenfunction expansion, and hence start by assuming that u has an expansion in terms of 


the eigenfunctions X (x). Thus, 


NITX 
u(x,t) = 2, u, (t) cos( (Eq. 2.10) 


eigenfunction expansion 


where the coefficients u(t) are functions of t. We now expand gq as 
q(x,t)= >  q,(t)X, (x) (Eq. 2.11) 


After substituting (Eq. 2.10) and (Eq. 2.11) into the 1-D inhomogeneous heat 
equation (Eq. 2.1), we obtain 


Xx 


Yu, (t)X,(x) =a, du, (2) |x (x) 1G (t)X,, (x). 


This yields the differential equation of the coefficients wu, (t), 


u, (t) =a, (2) hs (+74, (t) (Eq. 2.12) 


T 
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The initial condition for this equation is u,(t =0)=0. 
Now the problem is to find g,(t). We can use the property of orthogonality of 


eigenfunctions. By multiplying eigenfunction X,(x) which is cos(—") and then 


Xx 


integrating on both sides of Eq. 2.11, we get 


fala X, (adr = [Yg,OX,(x)-X, @dr= 4,(0 | X,2dae 


L, 1+ sours 


L Lc 
— ft ~ X= t)— 
an{)| —>——4 an()= 


L Ee Gey 


2 2 ' nITX 
Therefore, 4, (t) = L. q(x,t)X ,(x)dx = L. Le i i al 


To obtain qg,(t) 1s very expansive because the x—.x.(t) term makes the integration 


time-dependent. Here we use change of variables to make this calculation easier, 


x=x—-x(t) 


dx = dx 
thus, 
L,-X,(¢) _*3 re - 
q,(Q=— | he” cos(—(x+x,()))dx 
x —-X.(t) L, 


We introduce a constant d, which 1s the effective radius of the laser beam, and 


k ~2 
——_—x 
1) 


d is about 3 to 5 times of the standard deviation 6 of J,e . Instead of taking integral 


from x=-x.(t) to L.—x,(t) , using the assumption above to change the integral 


fromx=—d to x=d, since the laser beam is a very concentrated Gaussian distribution, 
so taking integral with respect to a large range is unnecessary; we only need to integrate 
within the effective laser range which is from —d to d. Figure 3 shows the effective 


radius of the laser beam. 
10 





—d () d x 


Figure 3. — Effective radius of laser beam, d as a multiple of standard deviation 6 


We also recall that the trigonometry formula of cosine in sum and difference 


formula, 
cos(u + v) = cos(u) cos(v) F sin(u) sin(v) 


then q(t) becomes 


(t)= = I ae vos) x eT asa 
qd, L 0 L L 


—d Xx xX 


—d x xX 


2(¢ a x (t) 
2% NX. NIX 
—-—| |JI.e”°  sin(—)dx |sin(——— 
: | [tc a | — 
The integral in the second term is zero since the integrand function is an odd 


function and the integration interval is |—d ; d|. Thus, g,(t) simplifies to 


x \ -d x x 


(t)= a: I aw cos) x et pladcaieud 
qd, L 0 L L 


Once q,(t) is found, we have to compute u(t) using the method of integrating 
factor. 


Recall the method of integrating factor in solving first-order ODE: 
1] 


yiQ)+a- yO)=r@) 
Here a is a constant and r(f) is a function of time. The first step of integrating 


factor is to multiply e on both sides of this equation: 


eH | y'(Q)+a- y(t) |= eH! .r(t) 


This equation 1s equivalent to 
d 

=a yO) = el ‘r(t) 
[ 


We can integrate on both sides of the equation and get 
t t 
e“ . y(t) = y(0) + | e“ .r(s) ds 
t=0 


at 


where s is a dummy variable. We multiply e “ to obtain 


t 
y(t) = y(O)- a ae | ae r(t) dt 
1=0 
Now, we can compute u(t) of Eq. 2.12 use the integrating factor method to get 


Ar 
K, 





u, (t) =a, -(=) u, (t)+2q,(t) 


























where the initial condition uv, (0) has been applied. 


. nx 

Once u(t) is calculated, we can find u(x,t) = DU, (t)X (x)= dU, (t) ae 

After we conduct several experiments in MATLAB, we find that the direct 
summation is very slow; the computational cost of this eigenfunction expansion method 
is extremely high in a 3-D nonhomogeneous heat equation problem. An alternative way 
to compute the series is to take advantage of the Fast Fourier Transform (FFT). In order 


to do so, we first change the cosine form to a combination of exponential forms. 


na na 
L L 
NIX e * +@ 
2 —s 


Xx 


Then we even extend u(x,t) from the domain (0,LZ,) to (-L,,0) and make it 


periodic with period 2Z, . Figure 4 shows how to expand the domain in FFT method. 





=” i 0 L 


One Period 


Figure 4. | Domain expansion: cosine expansion to FFT 
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So, we rewrite u as 


= hs + re 
u(x,t) = > u,(t) —+ 
n=0 
2 mx My (t) n > (0) 
= 'c,(t)e = where ¢,(t)=4 7 
to Aaa) ) , n<0 
2 


and then we apply FFT. We will provide more details on how to compute the solution 


with FFT in Chapter III. 
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Hil. NUMERICAL SOLUTION FOR A TRANSIENT, ONE- 
DIMENSIONAL TEMPERATURE DISTRIBUTION IN A FINITE 


ROD DUE TO A DITHERING LASER BEAM 


In this section, we review some numerical methods for solving the 1-D 


nonhomogeneous heat equation. 


A. THE CRANK-NICOLSON METHOD 


The Crank-Nicolson method is an implicit, unconditionally stable numerical 


method used to solve the heat equation. It is second-order accurate in time and second- 


order accurate in space. In 1947, Crank and Nicolson suggested an alternative way to 


utilize centered differences scheme [7]. To illustrate the Crank-Nicolson scheme, we 


start with centered difference and return to the Taylor series expansion for f(x, + Ax) and 


Ff (% — Ax) 


Ax) Ax) 
fore any=fay+ A pay BD preys! ) pO) +. 


a4 


Ax) Ax) 
Fa ax)= Fe) 2 prey BD prey AD prays. 


Subtracting (Eq. 3.1) from (Eq. 3.2) yields 
f (xt+ Ax) — f(x—Ax) = 2Ax f (x) += (Ax) fO(X)+... 


We can get f '(x) by reorganizing this equation 


_ f(xt+Axr)— f(x-Ax) 1 2: A) 
a ; (Ar) Vie Co eare 


=O(Ax)? 


f(x) 


This leads to the first-order centered difference approximation 


J (xt Ax) = f(x= Ax) 


f(x) * Te 


(Eq. 3.1) 


(Eq. 3.2) 


This approximation is with truncation error O(Ax)’ , an order of (Ax) and the 


approximation 1s accurate and consistent; the truncation error vanishes as Ax > 0. 
By adding Eq. 3.1 and Eq. 3.2, we have 
2 " 2 4 (4) 
f (x+Ax)+ f(x—Ax) =2 f(x) + (Ax) Ff (x) +7 (Ax) foo (x) +... 
Again, we can get the second-order derivative f"(x) by reorganizing this equation 
2 
_ f(xtAx)-2f(a)+ f(x—Ax)_ (Ar) 


(4) 
(Ax = fC (aX) +... 


O(Ax) 


"(x) 


and this leads to the second-order centered difference approximation 


J(xt Ax) 2 f(x) + f(x— Ax) 


f(x) ® (Ax? 


After we understand the second-order difference approximation, we can move on 
to the Crank-Nicolson scheme. The forward difference in time of temperature 


T =T (x,t) can be written as: 


or ~ T(t+At)-T(t) 
Ot At 


The Crank-Nicolson method may be interpreted as the centered difference 


A . tee, ofl T At. . ; 
around f+ = The error in approximating —( + = is O(At) [7]. Thus, we discretize 
t 


ee At. — 
the second derivative at oat a centered difference scheme. Since this involves 


functions evaluated at this in-between time, we take the average at tf and t+Ar. This 





, , _ OT °T 
yields the Crank-Nicolson scheme for 1-D homogeneous heat equation = = a, - ae 
t IX 
4 ages "7 y bls ee a IT* fe aes a OH ir ete 
i a ar i+] 5 Eg i=l ef i+] i ; + b= (Eq. 3.3) 
At 2 Ax Ax 
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We introduce the shorthand notation T(x,,t,)=T,, which is a matrix notation 


used to discretize x and t, because the Crank-Nicolson scheme involves six stencil 
points rather than a simpler stable method, and three of which are at the advanced time, 


as shown in Figure 5. We cannot directly march forward in time. Instead, we are left 
with a traditional system of equations of the form Ax =b , where X is the solution at the 
(k + 1) st time level and b depends on the solution at the k-th time level. Of course, it is 


not just any A, the matrix we have has a very special structure which 1s to be shown later. 


t+At 


x—Ax x x+Ax 


Figure 5. — Implicit Crank-Nicolson scheme [From 7]. 


Now consider a 1-D nonhomogeneous heat equation with a heat source f (x,t): 


0T OT 


Ot Ox’ 





+ f(x,t) ,O<x<l 


The BCs are insulated at both ends of the rod: 


or 
Ox 


_oT 


22) 26 
= ON 








x=l 


The IC assumes that there 1s no temperature rise with respect to the ambient 


temperature: 


T(x,0) =0 
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In this thesis, we study the case where the heat source 1s a dithering 1-D Gaussian beam: 


—(x-x,.(t))? 


2 
f (x,t) = Ihe ° HOen tos ye 





where J, is the intensity of the laser beam. As before, we discretize this 1-D rod into n 


parts, as depicted in Figure 6: 





Figure 6. _—Discretization of the 1-D rod 


Here in the Crank-Nicolson scheme and later in the matrix notation, the BCs are 
extended into ghost grid points x_, and x,,,. However, in this finite difference method, 
the BCs link temperature at x_, to the temperature at x, and the temperature at x,,, to the 


temperature at x. More specifically, the BCs are 








OT) Fst) ~The) _ 9 pk it 
ax Ax > 
T —T 
or _ wy Vow IW Pt) grt a7 for allk =1 
6) la Ax 


Applying the Crank-Nicolson scheme to the heat equation with a _ heat 


source f (x,t), we obtain 


TTS 7 a i ie —27" +T*, | Wb gran emer) Bik 


i+] 


l 
a2 a a eS Gri )+ FG rhea) 
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i=1,2,...n—I,n aad eee t, =kAt 


Here, 17; hare known, 7," are to be found. We start with a simple case i=1 and 


rearrange this equation: 
r aAt| T; -278+T. Ty 270 47%" At 
T;" = T* +—— Ss —— 5 LP oh + fst) 


where 7, =7\“ and T“"' =T,“"' because of the BCs. After changing some orders, we put 


T,*' on the left-hand side of the equation and 7," on the right-hand side; we get 


aAt aAt 


aAt Tk 
2Ax> 7 2Ax? ' 


2Ax" 


aAt 


yeas. = 
"DA? 


k+l k+l] __ ee 
Ty + Ty =t + 














T;" + LF ait eect te.) 


To simplify, we get 














c at ri" )- at (ae -(1- aAt r+ aAt 7) +L Ft) + Fst] ( Eq. 1) 


To put this in matrix form, we get 


at at qT 
I+ 20 9 Ax Pl oa 
2Ax AX Sixs I; 2x1 


aAt at (be At 
=|1- x} | +—| f(x,t,)+f(%.h, 
JAX? JAX L. kag 2 | f ( 1 i) f( Lek )] 


2 














Likewise, when i = 2 


T,' T,: T, — 27; +7 T,' - 3 eae Bi 
—— -9|R-nT oo + SLi Cent, \+ f (Xtra) 


After changing some orders, we put 7;“*'on the left-hand side of the equation and 


T,“ on the right-hand side 
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Ot ihe s{u-Se"|- Ot eae fo 
ON?! Ae 2Ar ° 
aN a) 4 | aM (Fi) 
aye {1-4 yen +S LF t+ f(s, t..1)| 
In matrix form, 
pes 
1 
- = (pace 7 ra ” ha 
2Ax Ax 2Ax 1x3 k+l 
T. 
3 3x1 
Bg 
1 
at a\t At At 
lar ae Jag |" /2 | tp erro) 
3 3x1 
When i=n-1 
TH _T* a T _9T* Pay fhe —27" + jie 
ee an <n +L, ity + f(%,_ 1? ti) 


Once again, we put 7;‘*'on the left-hand side of the equation and 7; on the right 


hand side 





AZ n-| IA? n 


Ay n-2 TEI = et Mis i 


(Eq. n-1) 


Ot 
apyeeie {1 ae are pot rani XN t, +f (x nt? t.1)| 


In matrix form 














‘ ; r To 
= aa | x| TH 
QAx Ax 2A" Ih kl 
T 
n 3x1 
yi 
ast . aAt aAt | | Ae 
_ At IT | +] f Out) +f Gratien) 
Fe Ax? 2Ax’ bd | ak 2 _ _ 
n 3x1 


Finally, when i=n 


Th Tk a =P ay. pel Ta 
At 2 


n+l ile nt a u 
—s re PSOne 


Here, T", =T* and T**' = | a because of the BCs. 


oa T, n+l 


Putting 7“'on the left hand side of the equation and 7;‘on the right-hand side 


leads to 














ON rest {ts at 1 |e 


wy n—l AY rhe saat | f(, t+ f Sst) (Eq. n) 


(eee 
2Ae "" (1 2Ae 


In matrix form, 
at aMt T* at at T* At 
— 1+ x| "I = |- x "| 4+—I[ f(x ,t,)+ f(x, 
| 2Ax° 2Ax* I. ad be 2Ax? I. i I 2 FC " r) IC " i) 


Putting all these n linear equations Fg. 1, Eg. 2,...Eg.n—1, Eg.n together in matrix 














form, we atrive at 


es 2 3 
AL nx = BL + i = Dax 
where 
jh (Be F(Qt d+ fst) 
hae ey F(X t+ F(t) 
=a ah — At 
— T —= — 
I 2 
Bi ‘tie F Garth +I Oa tea) 
re 7 T* 7 F(X t+ f(%G ste) - 
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— 0 0 
2A? —-2Ax 
wn teas § 0 
we Ae DAK 
0 4 
A=}, , 
0 
oN | oN aN 
0 ae D ee 2 
Ie Ae OM 
0 e- 10 a ee 
2A? AY I, 
Lo ee 0 
2X DAK 
0° eg 
B=| | 
= 0 
ONt ON aN 
a ae 
2NX Ae 2A 
0 : 4+ 2. i 
2X DIN Na 
i-s 5B + STS + EL Ph) + Pio 
ot aMt\ 4, Mt 4, At 
ae. tees (1-S)n Ae I, + 5 [ft +f Gta) 
._ . 





Ot at at 
FAY oeey as {1- Ae a Nf FA” T. +S UFC, fee Leer t,.1)| 


OMt _ x 





“{1- a Re + SUF )+ f(x Xn? 1) 


ZZ 


Note that A and B are triangular matrices. We can use the sparse matrix solver in 


MATLAB to create them and they are easily to be computed. In MATLAB, we simply 
use T =A\b to find the temperature rise with respect to time. 


Let us recall from Figure 2, the idea of a dithering laser beam heat source 


f (x,t) on a 1-D rod of Eq. 2.1 can be written differently: 














2 
=a, <4 +27 F(x,1), O<x<L, 
See 
(Eq. 3.4) 
I == %(0) : 
f(x) = —e , X(t) =x) +a-sin ,xX == 
2nd? period 2 


Where d is the standard deviation of the Gaussian heat source. Figure 7 demonstrates 
what the final temperature rise looks like for time = | using the Crank-Nicolson method 
with number of points n= 200. All the input parameters are listed in Table 1. The 
MATLAB code 1s attached in Appendix A. 





Table 1. | 1-D rod MATLAB input parameters for dithering laser beam 
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1-D Crank-Nicolson method with n=200 @ t=1 
1.1 





1.05 








Temperature 


- 
CO 
oO 





0.9 









































0.85 
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Figure 7. 1-D Crank-Nicolson method 


B. THE FAST FOURIER TRANSFORM (FFT) METHOD 


Recall that in (Eq. 2.10) we have shown that the solution of the 1-D 
nonhomogeneous heat equation with insulating boundary conditions can be expressed as 


eigenfunction expansion, where the eigenfunctions are cosine functions. 


Now we show how to use FFT to implement the eigenfunction expansion for the 


purpose of fast summation. This can be achieved in several steps. 


First, we even extend the initial condition from [0,1] to [1,2] and then make it 


periodic with period 2. After that, we map the interval[0,2] to [0,27]. By doing this, the 


Fourier series will only contain cosine terms. 
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Second, we apply the discrete Fourier transform to _ the heat 


e aig On a L)s 
2 Ot Ox’ —- 


Ou Cru 
|S = Fr) FL POs) 


where F denotes the Fourier transform operator. 


Now we need to calculate Fourier transformations of derivatives ofu(x,t). We 


begin by recalling the spatial Fourier transform of u(x,t): 
; { 2 | 
F | u(x, t) | = u(k,t) =— | u(x, tye”*dx 
2h 


Note that this 1s also a function of time; it is an ordinary Fourier transform with 
time ¢ fixed. To obtain a Fourier transform in space, we multiply e““ and integrate. 
Spatial Fourier transforms of time derivatives can be derived easily because the spatial 
Fourier transforms of a time derivative equals the time derivative of the Fourier 
transform: 

F ED = Jf ut) eee = AE Jacsnetnas = Ta 

Ot DT: OF O27 Ot 


For spatial Fourier transform of spatial derivatives, the method of integrating by 


parts can be used: 


imkx 


20 
F| _t | OU(%,1) jikee py. — UOC 
Ox QI Ox ee 





X= 2 Di 
—ikz E | u(x, near = —ikztu(k,t) 
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Here, we assume u > 0asx—-0 or 2z, and then the first term of this equation 


vanishes. 
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In general, the Fourier transform of the n” derivatives of a function with respect 
to x equals (-ika) times the Fourier transform of the function, assuming that 
u(x,t) > 0 sufficiently fast as x approaches periodic endpoint x > 0 and x > 2z. 


Likewise, Fourier transforms of a second derivative can be obtained: 


aed --iknF ee =(-ikz) u(k,t) =—k?nu(k,t) 
X 


Finally, we can conclude that "= (ik) Fat) =x at F(t) The 
t 
factor z here is due to the fact that we map the interval[0,2] to [0,27]. Using the 


method of integrating factor and multiplying e**" on this equation, we find 


2.29 du pe eas oe 
ea ke ute” F (kt) 
t 


We rewrite this ODE as 


dt 
Integrating both sides yields 


: (wel) =e! F(k,f) 


u(k,t)=u(k,0)e**' +e**' [ek F(k,s)ds =u(k,O)e**' + fe F(k,s)ds 
0 


0 
Here, we make the change of variables tr =t—s to conduct our following derivation. 


This equation becomes 


t 
u(k,t) =u(k,0)e*7' + | oe ek t—r)dt 

0 
In order to estimate this integral, we approximate flk,t —t)by a linear function 


(kt-1)= 7+ OE), 
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Then, it follows that 


u(k,t) =u(k,0)e "7! + | oe ek t—t)dr 
0 


© (ker + [ere Fuso ae fa 
0 


Z kn? iy l =k 7" £(k,0) — f(k,D) | —k? x? 
= kant ee: k°n°t 
=u(k, De" + (ks (1 e ———— oe *rdr 


A spe ”A 1 ee: 
= u(k,0)e* + p(k. (I-e k 


EET t eer, | (i-e**) 


t ken? koa" 





A ey ‘a 1 kt 
= u(k,O)e , + r(kst)aa(I-e : ] 


EASED 
[ 





1 l4+tk’n _», 
— —————_ @ 
k+r* kt nr* 

A.A, 


1-(L+0k? 2? jer 


k+n* 
Once u(k,t) is found, apply the inverse discrete Fourier transform to it to 


get u(x,t). 


In our MATALB code, as attached in Appendix B, tf is replaced by dt because in 
each step along the time direction we march dt. The first component in /, or 
h, corresponds to the case where k=0. We apply H’Lospital’s rule to find the value 


when k =0. More specifically, we have used the formulas: 


1 eay\ee* ,. 1-e*% (0 (ne 
fim Fi? (1 —e = lm ype] = lim —————- = t 








: 1-(l4tk?2’)e**" Pe Ie 0 ple ae" (0 
mo kg io ee eee oy) Dee 


Zi 


Figure 8 shows the result using FFT method in MATLAB where all the input 


parameters are from Table 1. 


1-D FFT method with n=200 @ t=1 
1.1 











Temperature 
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0.9 
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Figure 8. 1-D FFT method 
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c. COMSOL MULTIPHYSICS 4.0A 


COMSOL is one of the popular computer simulation software programs used to 
model and simply translate real-world physical laws into the real world in the virtual 
form using the finite element method. COMSOL is a commercial problem-solving tool 
that produces results quickly. However, it is more important to investigate the accuracy 
of the numerical results. We compare the solution differences between MATALB and 


COMSOL mn section D. 


Figure 9 is the result from COMSOL and there 1s no surprise that both results 
from MATLAB and COMSOL are very close. The step-by-step process of creating this 
1-D dithering laser problem using (Eq. 3.4) in COMSOL 1s attached in Appendix C. 


COMSOL 1-D 


Temperature 





Figure 9. = 1-D COMSOL 


a? 


D. COMPARISON ON A MODEL PROBLEM 


We compare the results from the 1-D Crank-Nicolson method and COMSOL and 
plot the relative difference in Figure 10, where the number of points is increased from 
200 to 2001. The relative error between MATLAB and COMSOL results is about 10%, 
which is really small and tolerable. Here, we only show the comparison between Crank- 
Nicolson method in MATLAB and COMSOL, and the comparison in FFT method and 
COMSOL are also pretty similar. After several tryouts, it is intuitive to see that the 
relative error goes down once we increase the number of points N, but it takes longer to 


do so. 


x 10° Differnce between C-N Matlab & COMSOL N=2001 


























Relative Temp difference 
Ce) 












































Figure 10. Relative error plot in 1-D Crank-Nicolson and COMSOL 


E. REAL PROBLEM SIMULATION (STEEL AISI 4340) 


After we have demonstrated that both MATLAB and COMSOL are giving us 
acceptable results, we can choose a specific material to simulate a real problem. The 
material we use in our 1-D plot comparison is Steel AISI 4340, a built-in material in 


COMSOL which has the material contents listed in Table 2: 
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Heat Capacity J/(kg*K) 


Thermal Conductivity ei W/(m*K) 


Thermal Diffusivity 


C, 
Melting Point i 1783 


Table 2. Thermal property of Steel AISI 4340 





The | laser we deploy has the inputs in Table 3: 


Magnitude of Gaussian source ae W/(m‘3) 
Effective radius of Gaussian Effective radius of Gaussian heat source source Eseseronoccominnernns [4 lowe [fm __ 


a 


Table 3.  1-D i laser input on Steel AISI 4340 





Figure 11 from MATLAB and Figure 12 from COMSOL depict the results of 
temperature rise using the contents from Table 2 and Table 3. Both results are very close. 
Figure 12 is from the Crank-Nicolson method and the FFT method yields a very similar 
result. It is clear that the maximum temperature rise decreases as the rotating period 
decreases; in other words, the higher the frequency is, the less the temperature rise 


increases. At t=ls, Figure 13 shows the maximum temperature rise of the steel AISI 
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4340 versus the frequency (reciprocal of the period) of the rotating Gaussian beam from 
1 Hztol100 Hz. It is observed that a higher frequency leads to lower the maximum 
temperature rise. It should be pointed out that our results are consistent with earlier 


analytical studies for the semi-infinite domain [5]. 


Figure 14 shows the temperature rise as a function of time at a fixed point 
x =0.75 of the material with dithering laser beam period=0.1 sec. The temperature 
increases as the dithering laser beam moves closer to the point and the temperature 
almost stays steady as the laser beam moves away. The overall line shape behaves as an 


increasing function of time. 
Figure 15 shows the maximum temperature rise versus frequency with two 


different heat sources/,. The higher J/,, the higher the temperature rise. 


1-D Matlab CN 
900 




















Temperature 
AN 
o) 
Oo 















































Figure 11. Temperature rise on the material of steel AISI 4340 using the 1-D Crank- 
Nicholson method 
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Figure 12. Temperature rise on the material of steel AISI 4340 using 1-D COMSOL 
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Figure 13. 1-D maximum temperature rise of steel AISI 4340 versus the frequency of 
the dithering laser beam 
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Fixed point x=0.75 temperature change 
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Figure 14. 1-D temperature change with time at a fixed point x=0.75 of steel AISI 
4340 with dithering laser beam period=0.1s 
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Figure 15. The maximum temperature rise of steel AISI 4340 versus frequency with 
two different heat source J, 
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IV. NUMERICAL SOLUTION FOR A TRANSIENT, TWO- 
DIMENSIONAL TEMPERATURE DISTRIBUTION IN A FINITE 
FILM DUE TO A ROTATING OR DITHERING LASER BEAM 


Now we extend our work in 1-D to 2-D. 
A. THE CRANK-NICOLSON METHOD 
We start with the 2-D nonhomogeneous heat equation 


Cu Ou 04 
eh eae +— f(x, y,t) , OSx<L, O<yK<L 
Xx y 
Ox” Oy ) K, 
~((x-x, (1) H(y-y, (1))”) 
2d’ 


X,y,t =—+e ’ 
IGN = 5 ap 





(Eq. 4.1) 





} XxX) = 
period 


L, 

2 

sf 2 L 
y.(t)=y, +b-sin 2 ,¥=— 
period 2 


X(t) =X, +cos| 





(Eq. 4.1) can be illustrated as in Figure 16. The heat source f(x, y,t) created by the laser 


beam is illustrated in Figure 17. 


Laser beam 


(0,L,) (L,,L,) 


counter-clockwise 


rotation 


(x,y, (0)) 


d 
(x(t), y¥.(4)) 





(0,0) (L,,0) 
Figure 16. 2-D schematic of the laser beam and the finite work piece 
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Laser Beam TEM 





2-D Gaussian Distribution 


—((x-x, (t))? +(y-y,. (1) ) 
2d? 





Heat source: f(x, y,t) =—"*“e 
f(x, y,t) rd? 


Figure 17. Laser beam creates heat source as a Gaussian distribution [After 2]. 


We discretize this 2-D rectangular domain as shown in Figure 18. In order to 


facilitate our presentation, we set L, = L, =1 and the thermal properties a, and K, equal 


1, thus the heat equation is simplified to the form = = Au+ f(x, y,t). 





Figure 18. Discretization of a 2-D rectangular region 
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The discretized point is denoted as (x,,y,), where i goes from | to n and 


j varies from 1 tom. Here, n and m are the total number of grid points in 


xand yrespectively. More explicitly, 


58 -(i-F Jax Mx= (L, ay. f=1 9 oh 


nN 





L, —0 
y, -(i-Jav, py {-9) j=12,..m 


The boundary conditions assume that the solid is insulated at the edges 
(boundaries) and the initial condition is that there is no temperature rise with respect to 


the ambient temperature initially. Thus, we have 








eu} _ au) _ 
Oils: OX 
BCs: 
eul _ du) _ 
OY los Oy}, 








ICs: u(x, y,0) =0 


We introduce the shorthand notation u;, =u(x,,y,,t,). Then the Crank-Nicolson 


2-D scheme of (Eq. 4.1) becomes 


k+l] k k k k k+l] k+l k+1 
U; ; —U; ; 7 ay Ui j —2u; TU: > Ui. j —2u; ; TU: 
At 2 Ax? Ax’ 
k k k k+l k+1 k+1 
% i Us a 2U; , +u; Us ia 2U; ; te 
2 Ay Ay? (Eq. 4.2) 
] 
+L POY t)+ LOY ote) | l<i<n, 1l<j<m 


k=0,1,2..... t, =kAt 


The BC can be satisfied by setting 
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Uy, = My, 
U, — Wrst 
y=, 
Um — ets 


: ; ; ie 
We can rearrange the temperature at the interior grid points to form a vector u as 


depicted in Figure 19: 





+k 
Figure 19. Putting the temperature of each point in a vector form wu after discretization 


We will convert the linear system (Eq. 4.2) into the following matrix equation: 
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>k-+1 —k 
where Aand Bare the matrices to be found, ft , are unknown and ft | are known. 


We rewrite (Eq. 4.2) as 








At 
k+l k k 
U; 7 es a JAX? Luh, — 2u; j TU; j 2Ay i, j+l 
At k+l k+1 k+1 At k+l k+1 k+l 
+ | i+, j —2u;,, +upi!, |+ 2Ay? I i, j+l —2U; +up" | 


+S TFG, Vit +f, jot) | 


Putting u‘"' on the left-hand side of the equation, u“ and f on the right-hand side, 





























we obtain: 
A A A A A 
: ak +{ 14 dl - : ur — a us a us — a Ted 
2Ax ia Ax” Ay’ |” 2Ax 4 QAy 7" DAy* 
At , At At | , At , At , At , 
— 2 AX? Wi, — Ax? — Avy? Je jv Ax Ui j DB Ay? it DAy? ei (Eq. 4.3) 


+S Osyth) + FOI pha] 


We start a simple case with i =1, 7 =1 first, 




















At k+l At At k+l At k+l At k+l At k+l 
“pat te ay? | aa Ot Day? Daye Oe 
At At , 


k 
+ Uy 9 
= 


Mat oa tot Daye Daye 








+ LF Wott fp Mota] 


e ee r 1 1 k 1 . e 
Using the boundary conditions (uj)! =, /', u;;' =u; and so on), we can simplify 


this equation 
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At At k+l At k+l At k+l 
lta DAy? | Ax? 2 DAy? 


-| At lu At , At , 


At 
“Bas? Bay? MO Dae Gaye tie tg LO Io) + FAL eh] 





k+l +k oo 
Recall the matrix form Au =Bu + f . We can obtain the corresponding 


components of A and B from this equation: 























A(1,1) =1+ 8 aa 
2Ax” 2Ay 
At 

AC, 2) =- 

U2) 2Ax? 

Miah 

2Ay 

B(L1) =1- At Ar 
2Ax” 2Ay 
At 

Bd,2)= 

U2) 2Ax? 

B(,n+l)=— 
2Ay 


A 
and we will put = F(X t+ FO, Votin)| into a component of a vector function 
called 7, 


When 7=2, j =1, (Eq. 4.3) becomes 

















iL. = = u u 
| vee yr 2 ~2,0 
2Ax 2Ay 2Ay “> 
=yk! 


At , aeae At , At , At , 





_At sia] At 4 en AP AP eg AP 








u,, +—— 4, , + — 4,5 + — 
a1 a a 2,0 
2Ax" 2Ay 2Ay* i 
4 


+ [Fes Yt I+ LO, Visti) 


To simplify, this is equivalent to 
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k+l At gat 


u = u = 
pm | 
*  QAx? 2Ay* 
—-—— pee 
A(2,n+2) 











| At k+1 





A(2,3) A(2,1) 
At At At At At 
U5 | ar TA th oF QAy? z, +S 1F Ca, Vist J+ F(%, ee) 


J 
B(2,n+2) 








B(2,1) 











u u 
| Dr al 
2AXx 








At 4 sek 








+ LF Yt )+ OG otea)| 


To simplify, this is equivalent to 
At At k+l At k+l At k+l 
Mai TA Uri — Ay cr 














A(3,n+3) 





A(3,2) 


k 


oa a 
3,1 2,1 a2 
2Ax" 2Ay" 
Re, some 








A(3,4) 
| At , 


B(3,2) B(3,n+3) 


When 7 =n, j =1, (Eq. 4.3) becomes 
At k+l At k+l] 
2Ay* os 


n, 
__ k+l 
Uy J 

















At , At a 


——Tu 
n,2 ,) ,0 
2Ay* 
ae 





At , At At 
u +|/1- Se 
Ay 








Ax? 


— 7 nad 
2Ax — 


=Uy | 


+L F A dote+ FO Mote) 


4] 


To simplify, this is equivalent to 











+ + u 
2Ax? 2Ay? |" 2Ax? 2Ay? " 





At At : MG gs A a 
[ | k+l uk yet 











A(n,n) A(n,n-l) A(n,n+n) 
At At |, At , At , At 
=||- — u., +———u Wet —| J (es Vick) es Vist 
| AY red nl Ay? nll 2Ay na Ty Lf ( we Viste d+ FXas Voth | 
B(n,n) B(n,n-l) B(n,n+n) 


Now i=1, j =2, (Eq. 4.3) becomes 




















y k+l 1 At At k+l At k+l At k+l At k+1 
- 7 Uy > a ee 7 Uo. — 7 U3 7 Uy 
Ay 2Ax" = 2Ay 2Ay 
Uy 9 


At , Coeaera Gere te A ke tal 
2Ay ~ 





1 Diep thy 


hy 3 


+S LF Vast )+ f(s Vasten)] 


To simplify, this is equivalent to 























_ At ie f l of. At ue At Ti _ At Th _ At us 
2Ar ~~ 2Ar Ay | Ayr Ayr 
A(n+1,n+2) A(n+1,ntl) A(n+1,(n+1)+n) A(n+1,(n+1)—n) 
At k At At i At k At k At 
a. gl {eK ==. f +—| f%™3Voot +S OG Votes 
i | 2A? rae day "Day On dah IAL, Yer) 
B(n+1,n+2) B(ntl,n+l) B(n+1,(n+1)+n) B(n+1,(n+1)-n) 


The same arguments hold for i=1,...,n and j =3,....m. 


For i=1, j =m, (Eq. 4.3) gives 





























At k+l At At k+l At k+l At | At k+l 
Dae em ta age |e age Hon Day? lent Daye 
At , At At | , At, At. Atk 
~ Dae om Ee Ay? [Mom aa loan * py? Mn Taye lo 


+S LF Vmot J+ I(%, Ynotea)| 
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To simplify, this becomes 


























At k+l At At k+l At k+l 
: D 2,m 77 2 [AL — 2 UW ml 
2Ax 2Ax” 2Ay 2Ay 
A(n(m-1)+1,n(m-1)+2) A(n(m-1)+1,n(m-1)+1) A(n(m-1)+1,n(m-1)+1-n) 
At t At At t At ke 
= D Uy, +| 1— 7 2 | 4m + 5 UW m-l 
2Ax 2Ax”  2Ay 2Ay 
— J 
ETE eee aol) B(n(m-1)+1,n(m-1)+1) B(n(m-1)+1,n(m-1)+1-n) 
At 
+S LEO, Vnoty) = eee Va 
For i=2, j=m,, (Eq. 4.3) takes the form 
At k+l l At At k+l At k+l At k+l At k+l 
2U3m TL ETI tg JY in ~ 2 im — 2 eo mt 2 Yo m1 
2Ax Ax” Ay 2Ax Ay —~— 2Ay 
=U m 


At , At , 











At , At At | , At , 
= a | ee 7 mn 2 49 m1 7 42 m1 
2Ax Ax” Ay 2Ax 2Ay’ —— 2Ay 
=U) m 
At 
+ 1F Co, Vnot d+ f(%; Vsti) 
To simplify, this becomes 
At k+1 l At At K+] At k+1 At k+l 
= ) 3,m TY] Lr 2 a3 2 2,m a 5 Uy in a 2 Uy m_| 
2Ax Ax? 2A 2Ax Ay 
A(n(m—1)+2,n(m—I)+3) A(n(m-1)+2,n(m-1)+2) A(n(m-1)+2,n(m-I)+1) A(n(m-1)+2,n(m-1)+2-n) 
At : At At |, At ; At : 
= 2 Us in an 2 > |4o.m + 2 Us in + 2 Uy nt 
2Ax Ax? 2Ay 2Ax 2Ay 
Ucx~—Y 
B(n(m-1)+2,n(m—I)+1) B(n(m-1)+2,n(m-1)+2-n) 


reais 
B(n(m-1)+2,n(m-1)+3) B(n(m-1)+2,n(m-1)+2) 


+ [fle Ymrb dt f(%3Ymote) | 


Finally, when i=n, j=m, (Eq. 4.3) becomes 
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We can put all the components of A and B together to obtain the two 


and B 


n-m,n-m) (n-m,n-m) ° 


matrices A 
+k —_ 
Unmx1 + f 


nmxnm 


k+1 


nmxl = B er | 


Uu 


nmxnm 


After obtaining A 


—k+1 


by u =A\b. 


is attached in Appendix D. 
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= Panwa , we solve it in MATLAB 


We use (Eq. 4.1) and the parameters in Table 4 to show the temperature rise of a 


model problem with the numerical points n =m =128 in Figure 20. The MATLAB code 





Table 4. = 2-D film MATLAB input parameters for rotating laser beam 
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(a) time= 10s, period= 1s (1b) time= 10.25s, period= 1s 
1 























time= 10.5s, period= 1s time= 10.75s, period= 1s 





























Figure 20. Snapshots of the temperature rise on a film of a model problem induced by 
a rotating Gaussian beam using the 2-D Crank-Nicolson method 


B. THE FAST FOURIER TRANSFORM (FFT) METHOD 


The structures of Fast Fourier Transform in 1-D and 2-D are similar; instead of 
using fft and inverse Fast Fourier Transform ifft commands in 1-D MATALB code, the 
2-D code uses fft2 and ifft2 to carry out the computation. Briefly speaking, one needs to 
even extend the problem from domain [0,1]x[0,1] to [0,2]x[0,2], make the problem 
periodic in both x and y direction with period 2, and then apply the Fourier transform and 
its inverse to obtain the numerical solution. Our MATLAB code 1s attached in Appendix 


I 
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Figure 21 uses the input parameters from Table 4 and returns a very similar result 
to Figure 20. The 2D FFT has a better performance than the 2-D Crank-Nicolson method. 
In other words, FFT can produce the result faster than the Crank-Nicolson method with 


the same number of numerical points. We will discuss the details in Section D. 


(a) time= 10s, a 1s (b) 


time= 10.25s, period= 1s 
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(Figure continued on next page.) 
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time= 10s, period= 1s 


time= 10.25s, period= 1s 





(C) time= 10.5s, period= 1s (d) time= 10.75s, period= 1s 





Figure 21. Snapshots of the temperature rise on a film of a model problem induced by 
a rotating Gaussian beam using 2-D FFT method in 3-D view 


C. COMSOL 


Figure 22 (a) is the result from COMSOL and Figure 22 (b) 1s the result using the 
Crank-Nicolson method in MATALB based on (Eq. 4.1) and Table 4 with the numerical 


points n=m=128 in xand_ y direction, respectively. The results are very close and we 


will do a point-by-point error analysis in Section D. The detailed process of creating this 


2-D COMOSL 1s attached in Appendix F. 
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Figure 22. Model problem comparison: (a) COMSOL (b) MATLAB FFT 2 method 


time= 1s, period= 1s 




















D. COMPARISON ON A MODEL PROBLEM 


As we have mentioned in Section B, FFT returns the solution faster than Crank- 
Nicolson; the efficiency comparison is shown in Figure 23. When N=2” points, we can 
see that FFT takes about 10° seconds but Crank-Nicolson takes more than 10° seconds to 
finish the computation. So FFT is about ten times faster than the Crank-Nicolson method 


for this test problem. 
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Efficiency plot: Crank-Nicolson V.S. FFT in 2-D code 
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Figure 23. Efficiency plot: Crank-Nicolson 2-D method versus FFT 2-D method 


In Figure 24, by using the result from Figure 22, we compare the relative 
difference in temperature rise with respect to time at a fixed point (x=0.5, y=0.5) in 
MATALB FFT 2 method and COMSOL using number of point N=256. The relative 
error is about 10° and this result is tolerable. The relative error goes down as we increase 


the number of numerical grid points N. 


x 10° Difference between MATALB FFT 2 & COMSOL N=256 





Relative difference 














| | 
O 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
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Figure 24. Relative error plot in 2-D FFT method and COMSOL 
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E. REAL PROBLEM SIMULATION (STEEL AISI 4340) 


We use Steel AISI 4340 as our test material with material properties from Table 2 


and the rotating laser we deploy has the following input in Table 5: 





Table 5.  2-D rotating laser input on Steel AISI 4340. 


Figure 25 depicts the temperature rise at different times within one period. It is 
observed that heat will not spread out quickly enough due to the properties of the Steel 
AISI 4340 material so the temperature at those points directly shined by the laser rise 
quickly. However, those points far away from the points hit by laser, for instance, the 


center point (x=0.5, y=0.5), has almost no temperature rise. 


51 


Lime =0,25 ponigd=i “s time= 0.5 period =1 





















































Hime=0.75 period =1 7c 
: 7 a 3417 


a0 

250 
1 00 

150 

100 
u| 

50 

= it 
_ La] 


Figure 25. Snapshots of the temperature rise on a Steel AISI 4340 film induced by a 
rotating Gaussian beam using COMOSL with period=1s 
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Figure 26 shows the temperature rise as a function of time at a fixed point 
(x=0.75, y=0.5) of Steel AISI 4340 with dithering laser beam period=0.1 sec. The 
temperature increases as the laser beam rotates close to the point and the temperature 
almost stays steady as the laser beam moves away. The overall line shape behaves as an 


increasing function of time, as expected. 
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Figure 26. 2-D temperature change with time at a fixed point of Steel AISI 4340 with 
rotating beam period=0.1s 


The quantitative relationship between the maximum temperature rise and the 
rotating frequency at time=Is is depicted in Figure 27; the maximum temperature rise is a 


decreasing function of the frequency (reciprocal of the rotating period). 


Max Temp rise VS Frequency at time=1s 











1 2 3 4 5 6 7 8 9 10 
frequency (Hz) 


Figure 27. 2-D maximum temperature rise of Steel AISI 4340 versus the frequency of 
the rotating laser beam 
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V. NUMERICAL SOLUTION FOR A TRANSIENT, THREE- 
DIMENSIONAL TEMPERATURE DISTRIBUTION IN A FINITE 
SOLID DUE TO A ROTATING OR DITHERING LASER BEAM 


We have already verified that COMSOL has returned a very accurate solution 
compared with other numerical methods in both 1-D and 2-D codes. Therefore, instead 
of writing a complicated MATLAB 3-D code, we use COMSOL to obtain the 3-D 
answer. Recall (Eq. 4.1) and impose insulated boundary conditions. We use Steel AISI 
4340 as our test material with material properties from Table 2 and the rotating laser we 


deploy has the following input from Table 6: 





baci 


Table 6. 3-D rotating laser input on Steel AISI 4340 


Figure 28 depicts the temperature rise at different times within one period. The 


hottest spot is where the point hit instantly by the laser as in the 1-D and 2-D cases. The 
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maximum temperature rise is 1468K. Those points far away from the area hit by the 


laser, for instance, the center point (x=0.5, y=0.5, z=1), have little temperature rise. 


Bre 27 ie, pecod= 1s Li fmpedsds, pert 15 


(a) 








(c) Hina =. Tis, pnd 1s La tire Tie, periode ls 





Figure 28. Snapshots of the temperature rise on a Steel AISI 4340 solid induced by a 
rotating Gaussian beam using COMSOL with period=1s 


Figure 29 shows the temperature rise at a fixed point (x=0.75, y=0.5, z=1) as a 
function of time. Figure 30 shows the maximum temperature rise of the whole domain as 
a function of time; the overall hottest spot is at (x=0.341, y=0.307, z=1) when time=0.64s 
with the rotating laser beam period=Is. Figure 31 and Figure 32 depict the temperature 
rise at fixed time=1s in different layers. It is observed that the heat does not spread 
downward quickly and there is almost no temperature rise 0.1m below the top surface 
shined by the laser. After several tryouts, materials have larger values in diffusivity than 
Steel AISI 4340 can make heat spread out faster and so make the temperature rise higher 


than Steel AISI 4340. 
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Figure 29. 3-D temperature change as a function of time at a fixed point of steel AISI 
4340 with rotating laser beam period=1s 
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Figure 30. Maximum temperature as a function of time with rotating laser beam period 
=Is 
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Figure 31. 3-D temperature rise at different layers from z=0 to z=1 
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Figure 32. Maximum temperature at different depth (a)z=1 top surface (b)z=0.99 
(c)z=0.95 (d)z=0.90 
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Figure 33 shows the maximum temperature rise at time=1ls when the period is 
reduced to 0.1s; the laser beam rotates 10 cycles. The maximum temperature rise 
decreases from 1754K of period=Is to 670K of period=0.1s. Therefore, the maximum 


temperature rise can be reduced by increasing the frequency of the rotating laser beam. 


Figure 34 shows the temperature rise at a fixed point (x=0.75, y=0.5, z=1) of 
period=0.1s as a function of time. The overall temperature rise oscillates, and its 
envelope behaves as an increasing function of time, but 1ts maximum temperature rise 1s 
smaller compared to the maximum temperature rise with period=Is, as illustrated in 
Figure 29. Figure 35 shows the maximum temperature rise of the whole domain as a 
function of time, and its envelope behaves as an increasing function as well. Figure 36 
depicts the quantitative relationship between the maximum temperature rise and the 
rotating frequency at time=ls. The results agree with earlier analytical studies for the 


semi-infinite domain [5]. 
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Figure 33. Snapshots of the temperature rise on a Steel AISI 4340 solid induced by a 
rotating Gaussian beam using COMSOL with period=0.1s 
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Figure 34. 3-D temperature change as a function of time at a fixed point of Steel AISI 
4340 with rotating laser beam period=0. 1 
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Figure 35. Maximum temperature as a function of time with rotating laser beam period 
=0.1s 
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3-D maximum temperature rise of Steel AISI 4340 versus the frequency of 
the rotating laser beam 
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VI. CONCLUSIONS AND FUTURE WORK 


In this thesis, both analytical and numerical solutions for describing the transient 
temperature rise induced by a moving laser in a finite domain have been developed. We 
have exploited several methods including the eigenfunction expansion, the Crank- 


Nicolson scheme, FFT and COMSOL. 


We have confirmed that the faster the laser rotates (1.e., the higher the frequency) 
the lower the temperature rise induced. In other words, to reduce the military’s 
vulnerability to high-energy laser weapons it is possible to let the object rotate or rock to 
minimize the temperature rise. The quantitative relationship between maximum 
temperature rise and laser beam rotation frequency can be used as a general guide for 
adjusting the speed of rotation of the object in order to prevent temperature rise from 


reaching the melting point. 


This thesis can be explored deeper in the future. Some future potentials 


endeavors include but not limited to: 


A. Increase or decrease the effective radius of the laser beam d in Figure 3 


and Figure 17 to analyze how temperature rise is affected. 


B. Increase or decrease the radii a and b of the rotating trajectory of the 
Gaussian beam in Tables 3, 5 and 6 to analyze how maximum temperature 


rise 1s affected. 


Cc Create 3-D MATALB codes and compare the results with COMSOL. 


D. Instead of using JEM, mode Gaussian distribution as the heat source 
illustrated in Figure 17, different transversal modes in a laser spot such as 
TEM,, or LEM,,, can be used to further seek the analytical and numerical 


solutions [8]. 


E. Thermal properties are assumed to be temperature dependent: this makes 


the nonhomogeneous heat equation nonlinear, but more realistic [3]. This 
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nonlinear model can be solved without mathematical background in 
nonlinear programming using appendix G, COMSOL code for 3-D 
simulation by imposing certain materials whose thermal properties are 


temperature-dependent, such as silver [9]. 


In COMSOL 3-D geometry, try cylindrical coordinate and spherical 


coordinate rather than Cartesian coordinate. 


Conduct an experiment and see if the theoretical modeling 1s accurate. 
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APPENDIX A. CRANK-NICOLSON CODE FOR 1-D SIMULATION 


N=256; 6 number of numerical points in [0, 1] 
dx=1/N; spatial step 

x={[TeN] ?=0.5) *dx; numerical grid 

u=zeros(N,1); solution at current time, u(j)=u(x(j)) 


oP o\0 


ole 


d=0.02; %& radius of Gaussian source 
a=1.0; 6 magnitude of Gaussian source 
uQ=zeros(N,1); % initial value 

u=u0; 

dt=dx/8; $ time step 

T=7 5/256; 

m=T/dt; %& number of time steps 


C=—(Orm|*dt? 


i) 


=n 7 ax 2: 

% forming the matrices 
G0=(0.5" onesinaz,1lj)2 Us5147> 
diL=0...0" Ones (N; 1) "xr; 

qd 1-045" 0Onesi(N, 1) “1; 


A=spdiagse ((=d_1, irdd0, -—d1l], [—-1,0,1),N-N)> 
B=spdiragse(([d_l, Lb=dd, d1i, [=L,0,1),N,N); 
VYO=0,570.2597Sin (L0*2* p10)? % location of Gaussian source at t*”{k=-1} 


£O>a*exn (= (e700). 27 (2* 02) ) sare (2 *%o1* a2); % Gaussian source at t*{k- 
ie 


= 


plot(x),u,'D=", linewidth’, 2.0) 

axe (0, 1,—0.04,0.16 |) 

drawnow 

tor k=1<¢m;, 
Vil=0.5Pr0s 257% Sin (10s Z* past (kel) %S location of Gaussian source at t%k 
fl=a*exp (-(x-yl) .*2/ (2*d%*2)) /sgrt (2*pi*d*2); %& Gaussian source at t%*k 


b=B*utdt* (£0+£1) /2; 6 right-hand-side of the linear eg 
tor wu’ tk} 

u=A\b; % solving the linear eq for u*{k} 

f£O0=f1; 


Sspause (0.1) 
plou(x,u, 'b=",; *linewidch*, 2.0) 
axis([0,1,-0.04,0.16]) 
drawnow 

end 
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APPENDIX B. FFT CODE FOR 1-D SIMULATION 


N=256; 6 number of numerical points in [0, 1] 
dx=1/N; spatial step 

x=[0O2N] **dx> numerical grid 

u=zeros(Nt+l1,1); solution at current time, u(j)=u(x(j)) 


oP ol 


ole 


d=0.02; S radius of Gaussian source 
a=1; 6 magnitude of Gaussian source 
uQ0=zeros(Nt1,1); 6 initial value 

u=u0; 

dt=dx/8; $ time step 

T=7 5/256; 

m=T/dt; % number of time steps 


C—=(Orm|*dt? 


o\? 


6 going to the coefficients of cosin expansion 

w=[u; u(N:-1:2) ]; 

Z=fft (w); 

cu=real (z(1:N+1))/N; 

= (OTN) oa) a Ze 

hi=[dt; (l-exp (=r (2:Nel) *dt))./r(23N+1).)> 

hnZ=(0.570E"2> (le (risen) *Orr 1) ..“ex0 (=r (2241) dt) ) J PZ INS) 2213 


= 


yO=0,570.2597Sin0 (L0* 2" p10); 6 location of Gaussian source at t*”{k=-1} 
£O>a*exo (=(e=70).° 27 (2 * 02) ) 7 sort (2 %01* a2); % Gaussian source at t*{k- 
i} 

w=[f£0; f£0(N:-1:2) ]; 

Z=f£ft (w) ; 

cf0=real (z(1:N+1))/N; 


S. 


DLO (x,y, t=", Litevio hn’ ,2.0) 
axis([0,1,-0.04,0.16]) 

drawnow 

ror k=Lim, 

Vil=O SPU. 2S Sint LO Ze piste (kek) )s % location of Gaussian source at t%k 
fl=a*exp (-(x-yl) .*2/ (2*d%2))/sgrt (2*pi*d*%2); %& Gaussian source at t%*k 
w=[f1; £1(N:-1:2) ]; 

Z=ff£t (w) ; 

cfl=real (z(1:N+1) )/N; 

update cu 

CuHCuUx *exp (=1r*dt) Fer, *nit+(er0-crl)7dc.=hz- 
erU=cr ls 

going back to the function 

Z=N* [cu; cu(N?:-1:2) |; 

w=1ifft (z); 

u=real (w(1:N+1)); 

spause (0.1) 

plottx,u, *r-*," Linewidch’® ,2.0) 

exis C0414 —O. 045.0. |) 

drawnow 

end 


oe 


o\? 
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APPENDIX C. COMSOL CODE FOR 1-D SIMULATION 


1. Open COMSOL 4.0a with 1D and hit = 
9 Model Wizard > [J Model Library 





select Space Dimension ao 


©3D 

(©) 2D axisymmetric 

© 2d 

() 1D axisymmetric 

@1D 
| OOD 
2. In Heat Transfer, select Heat Transfer in Solids (ht) and right-click Add selected then 
hit =. 

Add Physics 


> eS, Recenthy Used 
> * AC/DC 
> [8 Acoustics 
> Ses Chemical Species Transport 
> || J Electrochemistry 
> #26 Fluid Flow 
\\\) Heat Transfer 
| Heat Transfer in Solids (ht). 
BS Heat Transfer in Fluids ( a Add Selected 
24 Heat Transfer in Porows lwreore : 
« Heat Transfer with Surface-to-Surface Radiation (ht) 
ail Bioheat Transfer (ht) 
—e- Surface-to-Surface Radiation (rad) 
a Joule Heating (jh) 
> Go Plasma 
> Au Mathematics 




























! 


= 





+ 
Selected physics 





=a Heat Transfer (ht) 
3. In Preset Studies, select Time Dependent and hit Finish *. 


. fi Model Library — 





Select Study Type <= o> [ey] 
Studies 
@ 8) Preset Studies 


[= Stationary 


| i. Time Dependent 


> & Custom Studies 





Selected physics 


| = Heat Transfer (ht) | 
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4. Under Model Builder, right-click on Global Definitions and left-click to select 
Parameters, input Parameters as fae 















ydel Builder m=) \(1 Settings [ll] Model Library gC 
‘ ‘ | ‘Dm foc | | 
4 E Global Definitions | Parameters 
P| Parameters 7 Parumater 
a — Modell imag]! 
» & Definitions Narne . Vah | | 
JA, Geornetry1 d D.02{rn] 0.2m Gaussian Beam Radius 
# Materials | period i[s] is 
| — ht) , 125im| 15m nigel wicbon 
oi wf) 0.5(rn] 0.5m center of rotation 
‘ee WM 1 Hp sh nt tuce of ! 
M, sStep ls fhornO br peviod se i) L.O[Wi'm*3] Lyin magnitude of Gaussian Seam 
» [s, Solver Configurations 
— , spit : : ; . 
5. Right-click on Definitions and left click to select Variables, input Variables as 
following: 
del Builder 7 = 5 || ei Settings & Ql Model Library 
1-D.mph (root) eee 
= Global Definitions » Variables 
= Model (mcd) Geometric Scope 
= Definitions 


a= Variables 1 Geometric entity bevel: | Entine model 
LL View 1 


aN Geometry 1 + Variables 
€ Materials 
Wy Heat Transfer (ht! Name | Expression Unit 
3 Mesh 1 xe w)+a*sin(2*pi"t/period) m 
ea Suny] f D/sqrt(2"pitd* 2)[L/m]*exp(-(x-nc)*2/(2%d"2)) Wim? 


6. Right- click on Geometry and left-click to add Interval. The Left endpoint is 0 and 





the Right endpoint is 1. Left- click the build all button 

‘Ty Model Builder’. ~ = ©){#%2 Settings» (0 Model Library 

4 (9 1-D.mph (root) 
» E Global Definitions 
4 — Modell (mod!) + Interval 
= Definitions 

a P\ Geometry 1 Nurnber of intervals: | One ’ 
iA Interval (i) Left endpoint: 0 m 


wf Form Union fir 


to create SADE 








a & Materials Right endpoint: 1 m 
We Right-click ‘on Materials and left-click on add Materials. Make all the values in 


Material Contents equal 1. Right-click on the “line” in Graphics and left-click on it to 
select. Make sure that “ 1 ” is under the selection. 
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em Hes 
B&800 875 


Eude 7 E vl ite IB rredalitury 
ij Pega ea Pe 

Eie-al Te iitars * Materist 

Mocel L meal! Seninetr ke Eine 

B Cdwdcel 


fs Genres | fdeonibk ey ect (Caran 


fe Intereal 1 [2 ; : 
af Fore Lipebary peor [Para 2 
= meus 
OF ie Dafoe 
Hest Teeite Fc 
a Met] 
dads 
ih Sep Lee a) 
a, Toles Doar qunati, 
j faut 
© On Sas ries reais 
[oF Cerise Ve aer Eau: fe eeran 
=f Tarim Berwocheret, 
fi. ID Pat reve hecwrregwerhsdee 
Rapa Said ache 
© Pratl Pista edt rt br Wel ha, 
OS sen! Gar bebcacda tr 
cS Pert 





___ fropacty lies | the 
a Thaveal epadurialy © 
a Teac chit l 


of Fertcupacty woderi pu ts L 


8. Right-click on Heat Transfer and left-click to add Heat Source. 


1 Od @ | 8 el 


A) el Saphect 


c cl be + Be id gE Ce wT ag Vu 1 
SS Mew 25 8 Prego | Renuka ” 
| Ware _| depart [to Pentre ent 
We 6 ime 
kya ck 
Niet] dex 


For General 


source Q, Select User defined and Put “f’, which is the heat source defined from the 


Variable. 
a 0S 1-D.mph [reot} 
+ = Global Definitions 
a — Modell (mod) 
é ay Geome#try 1 
FA Interval 1 (ul) 
Form Union (fin) 
a #8 Matenals 
> Self Defined 
a GM Heat Transfer (ht) 
& Heat Transfer in Solids 1 
% |) Thermal Insulation 1 
= Initial Valores 1 
) Hest Source 1 
a 3 Mesh1 
pa Size 
>» [F Edged 
a & Study 
ba, > Step 1: fron 0 te period sec 
b fre. Solver Configurations 
a Results 


' Heat Source 


Domains 


Selection: Manual ¥ 


i 


pt 
x 1+ 


* Heat Source 


@ General source 





f Win 
-) Linear souree 


Q=q.T 





9. In Initials Values, make Temperature equal 0. 
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del Builder = C1) §i\Settings ». [ll Model Library SS 


1-D.mph [root! ™ Initi 
4 Initial Values 
= Global Definitions 


e Model 1 (mad) Domains 
+ = Definitions 
F A, Geometry 1 Selection: | All domains 
\ Interval 1 (it) ; } 
Form Union jfin! 
a % Materials 
» &8 Self Defined 
él Heat Transfer /hi! 
=~ Heat Transfer in Solids 1 
=) Thermal Insulation 1 
Initial Values 1 
-~ Heat Source] 7 Initial Values 
a G3 Mesh1 
Al Size 
b (2) Edgell T OK] K 
oa study] Surface radiosity: 
[d, > Step 1: from 0 to penod sec 
» Pes Solver Confiqurations Jo W 
10. Right-click on Mesh and select Edge, in Edge under Element Size, select Custom and 


Ellto build 





I 
[ 1 an” 





Temperature: 





make Maximum element size equal 1/200. Then select the build all button 











the mesh. 
mene — 7 = © || #44) Settings s. [I] Model Library’ qja> 
4 © 1-D.mph (root) " : 
b Global Definitions as vIZe 
4 — Madel1 (modi) Geometric Scope 
= = Definitions 
F ray Geometry 1 Geometric entity level: | Entire geometry 
Ay Interval 1 (i) 
Form Union (Fin) Alement Size 
a & Materials 
> &8 Self Defined (© Predefined | Extremely fine 
| = | Heat Transfer (ht) @ Custom 
S| Heat Transfer in Solids 1 - 
4 | Thermal Insulation 1 * Element Size Parameters 
Si  Initial Values 1 
) Heat Source 1 Maximum element size: 
a & Mesh1 1/200 - 
Ad Size - 
4 fer Edge1 ["] Maximum element growth rate: 
Ag Size1 14 
4 Ge) Study! [|] Resolution of narrow regions: 


i. > Step 1: from 0 to period sec | 
}1 


11. | Under Study in Step, select Range button ' , under Entry method, select 
Number of values, start from 0 and stop at “period” which 1s defined in the parameters. 


Put Number of values to be 20. 
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34 Settings ~_ [J] Model Library 





1-D.mph (root) ‘ 
= Global Definitions i lime Dependent 
— Model 1 (mad) + Study Setti 

= Definitions a am ae ear 

iA, Geometry 1 Times: range(0, (period[1/s]-0)/19, period[1/s]} 5 


A Interval 1 (ct) 
Form Union (fin) 








2 Materials 
Entry method: 
& Self Defined ntry | 
=a Heat Transfer (ht) Gant 
es |. a se 
= Heat Transfer in Solids 1 ae eye 


S) Thermal Insulation 1 
=) Initial Values 1 Number of values: 20 


| 
\ Heat Source 1 Function to apply to all values: 











Ea Mesh1 
= |_Replace | 
R 
[=] Edge - 
fax Size 1 a - —_ 
€3 Study 1 


[,. > Step 1: from 0 to period sec 
Tre. Solver Configurations 
Fa Solver 1 


Saf Compile Equations: from 0 to 


ui” Dependent Variables 1 
I, Time-Dependent Solver 1 Mesh 





4) Direct 
‘4 Advanced ~ Physics Selection 
a= Fully Coupled 1 Physics interfaces: 
Direct 1 Heat transfer (ht) 
Results _ 
=: Data Sets 
* Derived Values 
Tables 





1D Plot Group1 


Report Use in this study 
&) Player 1 Discretization: [Physicssettings § ——-—s—r—_——(sO | 


12. Under Time-Dependent Solver, make Relative tolerance to be 1.0e-6. Left Click 








= 
on —— to compute. 
corded Fhuilchier <I TYam)] em (“oO 
@ 1-D.mph (rocet} . 
Sabon Ae x, Time-Dependent Solver 
— tdodel 1 fot i) ~ General 
= Definithons 
. Geornetry 1 Tiers: ranged, (perked | 1s] 01/159, period[ 1s) = | heel 
Po, Interval 1 (ci) 
cma Relative tolerance: 1 Oe-f | 
= eames 4 + Absolute Tolerance 
| Detinec 
S% Hest Transfer (he) - Time Stepping 


S) Heat Transfer in Sealids 1 
5”) Thermal Insulation 1 
S Initial Values 1 b Output 
a Heal Sasuree 7 
Sad Meshi 
“4 Size 
Fiore 1 ‘ se 
24 Size 
fr. » Stepr Ds freer O fer perigee sec 
.. Sober Comiguratyons 
Fé) Solver d 
aT Ceempile Equestioees: fier 0 the 
wae Depordont Variables 1 
[¢, Time-Dependent Sohrer 1 
[a] [Threct 
“hy Aucthecel 
mae Fully Coupled 1 
Direct 1 


13. Finally, under Results, select Line Graph under 1D Plot Group and select time to be 


F Results While Solbvineg 


B Adsricendl 
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1 only, we should be able to see the result like Figure 9. 











4, Mieabrl Fachin iat = een, Liweilema, | GL | ll toate, aacRmiwe|m 
gga =. Line Graph : COMED 0 
2 tes! fetes nca | 
— Mek | oe [Linke 
© [Debwerry 
err; | fale ee] — tm 
bere I 
r ot Lie = ‘erie ee 
© bite us fas : tie} 
edad tl oo i 
Le Tracked ht LTT not 
| oe ie re 
A) Thered meketics J Geil 
WP) ieee eee: eee 4 
FP Kira Seer | ATI | 
Si kis Fb ae L is 
nr I Fin 
ES bebe See re + 
eer! 
Bhan I ee ee | i” 
Siig: ire Ode: i . 
"hy Seid Coe aor : ond. 
Fi, Soherl = 
og Corey Copier: Peo Be 
on Toe bea * tua 
jo Pre: Deperdret ahve | J 
3 eer nl 
bully tmgpireil 4 lat Tera e-. & : ; 
tent t | ie ia 3 is ms iu oJ Li oa i 
ly Be ren 1 
gj Cae ime fi 
fo Ci ell Pulte iia Vira = Tan FicrtT 
Tek ie CT 1 iT 
LO Fe Gg J L fais ead ie. 
», Log  Daectiptiors Coreebete sued ok 0 eee 
eat eee ery 
Pheer] = 
ff seers! 
EB rest nm PL ae 
Phra 


14. Further data analysis can be done under Report, such as generate a movie and export 
data and further compare results with MATLAB. 
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APPENDIX D. CRANK-NICOLSON CODE FOR 2-D SIMULATION 


N=128; %© number of 
M=128; S number of 


x= 1eN)] 0.5) *dx; 
y—( (LEM) “-O..5) ~dy; 


numerical points in x direction 
numerical points in y direction 


ie) 


6 Spatial step 


fe) 


6 numerical grid 


[xx, yy]=meshgrid(x,y); 


u=zeros (M*N,1); 
d=0.02; 

a=1; 

u0Q=zeros (M*N,1); 
u=u0; 


= 


dt=0 .5* dx; 
m=11; 
t=(Osm) *dt; 


=— 


Pi=ac7 de 2: 
r2=dt/dy%*2; 


= sparse (ll |y ily tl] 
=. Sperse (ly Li -i] 


is 
=1; 
(k,k)=1l4+r1/2+r2/2; 
(k,k+1)=-r1/2; 

(k, k+N) =-r2/2; 
(k,k)=1-r1/2-r2/2; 
(k, krl)=ri/2; 

(k, k+N) =r2/2; 


=N; 
(k,k)=14+r1/2+r2/2; 
(k,k-1)=-r1/2; 
(k,k+N)=-r2/2; 
(k,k)=1-r1/2-r2/2; 
(kh k=l) =f1/ 2: 
(k,k+N)=r2/2; 


FOr 2=Z2: (N=L), 


fe) 


Ss SOlUciCOn at curren: Eame 
radiue Of Gaussian scurce 
magnitude of Gaussian source 


fe) 


— I1nttial vyelue 


ole 


o\? 


fe) 


% time step 
% number of time steps 


/N*M, N*M) ; 
/N*M, N*M) ; 


% forming the matrices 
he bottom boundary 


A(i,1i)=l+rl+r2/2; 


tei—ljye=r 1/2? 
i,i+tl)=-rl1/2; 
i,i+N)=-r2/2; 


ty ial) =21/2> 
1p t+ Ly=Lil7Z? 


A ( 
A ( 
a 
Bi(aj,4 j=lerl-22/2; 
B ( 
B.( 
B ( 


1,i+N)=r2/2; 
end 
6 the middle layers 


75 


Lor 


end 


fe) 


Lor 


end 


J=2:(M-1), 
for 1=2:(N-1L), 
ka (=) Ne 
A(k,k)=l1+rl1l+r2; 
A(k, k-1)=-r1/2; 
k, k+1)=-r1/2; 
k, kK+N) =-r2/2; 
k, k=N)=-712/2:; 


end 


% left boundary 


J=2* (M=L), 

k=1+ (4-1) *N; 
A(k,k)=1+r1/2+r2; 
A(k,k+1)=-r1/2; 
A(k, K+N) =-r2/2; 
Atk, k=N)=—r2/2> 
B(k,k)=1-r1/2-r2; 
B(k,k+1)=r1/2; 
B(k, kK+N) =r2/2; 
B(k, k-N)=r2/2:; 


6 the right boundary 


ror 


end 


fe) 


J=2: (M=1), 
k=N+ (j-1) *N; 
A(k,k)=1+r1/2+r2; 
A(k, k-1)=-r1/2; 
A(k, K+N) =-r2/2; 
A(k, K-N) =-r2/2; 
B(k, k) =1-r1/2-r2; 
B(k,k=1)=ri1/2; 
B(k, kK+N) =r2/2; 
Bik, k=N)=f2/ 2; 


* the upper boundary 


k=1+N* (M-1) ; 


A(k, 
A(k, 
(k,k-N)=-r2/2; 
(k,k)=1-r1/2-r2/2; 
( 

( 


k)=1l4+r1/24+r2/2; 
k+1)=-r1/2; 


k, k+1)=r1/2; 
k, K-N) =r2/2; 


N 

k 

kal) Sari/7 2 
k=N)==712/2; 
k 
k 
k 
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EOr 2=—22¢(N=L), 

k=i+N* (M-1); 

A(k,k)=1+r1+r2/2; 

A(k, k-1)=-r1/2; 
A(k,k+1)=-r1/2; 
A(k, kK-N) =-r2/2; 
B(k,k)=1-rl1-r2/2; 
B ( 
B ( 
B ( 


end 

6% build the vector f£ 

f=zeros (N*M,m); 

fl=zeros (N,M); 

for k=l1:m, 
3CH—U.01r 0s Zo COSt LU 2% pack) )s 
Vc=0.570.29*5S1in(10*2Z* pit (Kk) ); 
PE=e exo (= (xe =xC) 27 (2d 2) Sno) a 27 2a) 7 (2p 2); 
£(:,k)=reshape (ff,N*M,1); 

end 

figure (2) 

drawnow 

ror k=l: (m=1),; 
b=B*tutdt* (£(:,k)+£(:,kK+1))/2; 
U=A \b- 
pause (0.2) 
temp=reshape(u,M,N); 
h=surf (xx, yy, temp) 
set (h, 'edgecolor', 'none', 'facecolor', 'interp'); 
SView (2) 
drawnow 

end 
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APPENDIX E. FFT CODE FOR 2-D SIMULATION 


clear 

N=128; 6 number of numerical points in [0, 1] 
dx=1/N; % spatial step 

x= [O2N]' *dxs 6 numerical grid 

[xa, ya]=meshgrid(x,x); 

u=zeros (N+1,N+4+1); 6 solution at current time, u(j)=u(x(j)) 
d=0.02; %& radius of Gaussian source 

a=1; 6 magnitude of Gaussian source 
u0=zeros (N+1,N+4+1); 6 initial value 

u=u0; 

dt=dx/8; $ time step 

T=12/128; 

m=T/dt; %& number of time steps 


t=[O2m)|*dt; 


oe 


oe 


going to the coefficients of cosin expansion 
[ 


w=[u; u(N:-1:2,:)]; 
w=[w,w(:,N:-1:2) ]; 
Z=f£fFt2 (w) ; 

cu=real (z(1:N+1,1:N+1))/N%*2; 
Na=[0:N]; 


[Nx, Ny]=meshgrid(Na, Na); 

LS (NX OL) a 2 NYy* pi) «ZF 

r(1,1)=1; 

hi=( 6x0 (=r*dt))<7/r- 

h1(1,1)=dt; 

ho (la(ee ort) .*exo(—=2"dr))./ fs 2 
NZ. (gd. H040" C0 2s 

r(1,1)=0; 


= 


[e} 


xU=0,570,250"COS (L0*2* p10) 3 


yYU=0.570..259*sin(LO*2*pa*0); os LOCaLLOn Of Gaussian source at Lt”, k=1} 
fO=a*exp (—((xa-x0) .*2+(ya-y0) .*%2)/(2*d%2))/ (2*pir*d’2) ; = Gaussian 


Source at ct k—1} 

w=(TOe TOC f—LsZ, 2). ]3 
w=|wy,wlty,N?t=—LtiZ) ]> 

Z=f£Ft2 (w) ; 

cf0=real (z(1:N+1,1:N+1))/N%2; 


fo) 
0 


surf(xa, ya, u, 'edgecolor', 'none','facecolor', 'interp') 
axis ( (0, Ly OU, 2, =—0U.04,0.4)) 
Caxis([0, 0.41) 
view (3) 
drawnow 
ror k=L?m, 
x10, 570.20" COs (L022 *pistec.(krh)); 
ViHO0. 570.2507 s1n (10424 pie (keh) 7 %* location of Gaussian source at t%k 
fl=a*exp (-((xa-x1) .*%2+(ya-yl) .*%2)/(2*d%2))/ (2*pir*d%2) ; % Gaussian 
source at t%*k 
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w=(tly £LLi(N¢HlLi2,:) 13 
w=[w,w(:,N:-1:2) ]; 
Z=f£FEt2 (w) ; 
cfl=real (z(1:N+1,1:N+1))/N%*2; 
update cu 
CU=Cl.* exp (=1*dl) +Crl “hike cr0=crl) 7c. he; 
cf0=cfl1; 
k 
6 going back to the function 
Z=N°2* leus. CUINt 122,24). 13 
Z=| 77. oy, Ni=lL22) )\s 
w=1fft2(z); 
u=real(w(1l:Nt+1,1:N+1))/; 
Sspause (0.2) 
surf(xa, ya, u, ‘edgecolor', 'none', 'facecolor', 'interp') 
exie(t0, ly UO, L =0.04,0.4]1) 
caxis([0, 0.4]) 
view (2) 
drawnow 
end 


oe 


80 


APPENDIX F. COMSOL CODE FOR 2-D SIMULATION 


2. Open COMSOL 4.0a with 2D and hit 
V3) Model Wizard > . o Model Library’ “8 Material Browser =O 


elect Space Dimension So 


Lit 
2 


) 2D axisymmetric 


) 1D axisymmetric 


(ee eo ty t 7 
Co 


= 
co 


2. In Heat Transfer, select Heat Transfer in Solids (ht) and right-click Add selected then 
hit =. 
Add Physics 


> eS, Recenthy Used 
> & AC/DC 
> ©] Acoustics 
cb Zee Chemical Species Transport 
> |] Electrochemistry 
> #22 Fluid Flow 
a \)\) Heat Transfer 
Heat Transfer in Solids (ht) | 

= Heat Transfer in Fluids ( .* Add Selected 
Heat Transfer in Porous hrreere 
r Heat Transfer with Sire tes Seats Radiation (ht) 
ail Bioheat Transfer (ht) 

<e- Surface-to-Surface Radiation (rad) 

ea Joule Heating (jh) 
> Go Plasma 
> Au Mathematics 








Selected physics 





= Heat Transfer (ht) 


3. In a ues select Time Peso and hit Finish A 


i tae Model Library’ O 





>Select Study [ype = o> || 
Studies 


@ 8) Preset Studies 

[= Stationary 

[P, Time Dependent | 
b> &) Custom Studies 











Selected physics 


| = Heat Transfer (ht) | 
8] 


4. Under Model Builder, right-click on Global Definitions and left-click to select 
Parameters, input Parameters as following: 

































"Ty Model Builder. | o #4 Settings », [ll] Model Library| & Material Browser| 
9 2-D.mph (t ‘ | 
' E he | STEUSUT 
“| Parameters |» Parameters 
() Modell mod) 
ee Study 
Results 0.02 m radius of Gaussian Source 
penod 15] 1s rotating period of heat source 
a ().29[m] 25 m racius of rotating trasectory in x direction 
b 0.25[m] 025m ragius of rotating trayectory in y direcban 
vil} 0.5[m] 0.5m center of rotation in x direction 
yi 0.5[m] 0.4m center of rotation in y direction 
i 1Des[Wim"3] LOEB Wim? Intensity of heat source 


5. Right-click on Definitions and left-click to select Variables input Variables as 
following: 


Gr fel Bui ler, 7 = Gi] 43) Settings odel Library Material Browser 
"TT Model Builder”. 7 “ Gi || es [i] Model Library 8 Material Browser | 
¥@ 2-D.mph (root) 

= Global Definitions 


()) Model 1 (mod]) Geometric Scope 
= Definitions 


a= Vanables 1 Geometric entity level: Entire model 


(& Boundary Syster 





a= Variables 








»L4 View 1 + Variables 
'F\ Geometry 1 
Materials _Name Expression Unit 
&@ Hest Transfer (At) xc +a" cos(2*pi*t/ period) mn 
E Mesh 1 yc y+ b*sin(2*pi*t/ period) m 
fea) Study 1 q TOs (2*pitd* 201 /m*2])*exp(-([x-xc)*2+(y-ye)*2)/(2"d*2)) Win! 
Results 


6. Right-click on Geometry and left-click square and make side length 1 m. Click 























build all button ‘«! to create a square. 
AED Mochel Builder > moo [ He Settings [DD Mosel Liteary | 2 Material Browser | awml/aua 
S 2-D.mph (reat! a 
E Global Thetis. Oo ela 
Model t imoxty) = Object Type 
& Definition: 
, Geornetry 1 Tree | Sclid - | 
Co] Squared (sq) 
af] Form Union ifin|| ™ Size 
M2 Materials Es 
SM Heat Transfer fhe) ieee 7 
ue bilesh 1 | |= Pest pain 
eh Sthincky 1 ; 
Results Bac Conver =| 
Pa a . mi 
i: oO ITI 


Rotation: oO dod 
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T Right-click on Materials and left-click on add Materials, add Steel AISI 4340. 
Right-click on the “square” in Graphics and left-click on it to select. Make sure that “1” 
is under the selection. 
‘TT Model Builder 
4 1 2-D.mph (root) 
> (& Global Definitions 
4 [) Model 1 (mod]) 
b> = Definitions 
» YA Geometry 1 
4 @ Material- 








— qi Model Lit 
® Materials 





Wak ari DY UL Stage bec Laberge res ao} » Uo e-|S S544 @ ee eo 
8 2Oo7p ral - j * | 
? ora Derdew ly 
Deke brent Seeciic Sooper 
we Defadar. 7 
 Cacrery 1 (pomeraches: clay bevet Lee + 
* lai maa: a 
Gel Teta L  # 
OB beet Saecter Pa fr - 
Masi - 
mn thay 
iy Resube 
dbase cacti 
Se Pqube 
Thetis heresy 
Chesorageabe Fd = 
ot betore 
Pheer niche, 
it. Ketel 
+ 
1 I I I I I 1 I 1 
Miatroad Comberety a a a a a eT ae: bo | / 
[oh Si pcs ee peewee Cy isleyg el | bvig*t in : wae Said bo Sh a AS 
| ft Senaity Pm ee a iain ; 
ot Ihe! cor ducretr k Afi] Wire 


8. Right-click on Heat Transfer and left-click to add Heat Source. For General 
source Q, Select User defined and Put “q”, which is the heat source defined from the 
Variable. 
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1 ) ils ; = = Fi] see Settings a Model Litsrary | SP Material Brewer | Gl 








# © 2-D.mph (root) rH ' 
: | eat Source 
e EE Global Definitions , 
> € Definitions. 
> FA Geometry 1 Selection: | hlanual ~| 
> GP Materials - 
3 . | 1 ™ + 
a = | Heal Transher (hi) 
& Heat Transfer in nH = 
& ) Thermal Insulat 
2 Initial Values 1 
)* Heat Source 1 
b 2 Shudyl 
ct Ge Results ~ Heat Source 
® General source 
O | User defined ~| 
qj Wim 


© Linear source 


9. In Initials Values, make Temperature equal 0. We assume there is no temperature rise 
in the beginning. 








"TD Model Builder ~~ 54) Settings) Ml Model Library | S$ Material Browser g@-a 
a © 2-D.mph (root) = Irniti 
= Initial Values 
. ( Global Definitions ~ 
af) Model 1 (mod) Nema 
» = Definitions 
" A Geometry 1 Selection: | All domains - 
>» §B Materials : 7 
4 Sm Heat Transfer (ht) sa 
&") Heat Transfer in ln = 
S| Thermal Insulati x 
= Initial Values 1 
Heat Source 1 
> G3 Mesh 1 
> BS Study 1 
Results + Initial Values 
Temperature: 
T OK] K 
Surface radiosity: 
Jo Wim 


10. Right-click on Mesh and select Mapped. Right-click on Mapped and select Edge 
Groups. Please select domain and add each edge group. Under size, select Custom and 
make Maximum element size equal 1/64 or any 1/2“N where N 1s an integer to compare 


: to build the mesh. 





with FFT 2 method in MATLAB. Then select build all button 
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TH Edge hroaups 


SO Thecrresal Ireruluti 
B™) Initial Vabucs 1 




















al Size 


i ean Sy = : 


— 


Geometric eniity loves |2ntne geametry 





©) Predefined [ Homa = 
Cueto 

© Homectit Size Parameters 
Vl Macmum «eect sor 
15a om 
7] Minimee element are 


RnR im 
] Mac mum @ eect qiach rate: 

q 
42 | 
7) Reechinien of runcature: Wtf - 
ke 





7] Reeclinicn cf nanices rep: 


ee al 


11. | Under Study in Step, select Range button’ . Under Entry method, select 
Number of values, start from 0 and stop at “period” which 1s defined in the parameters. 
Put Number of values to be 31. 
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"IT Model Builder F448 Settings \. fl] Mode! Library| % Material Browser] ZL = G 
V9 2-D.mph (root) 


= Global Definitions i, Time Dependent 








P| Parameters ~ Study Settings 
{)) Model 1 (mod) oe 
= Definitions Times: _range(0,(period[1/s]-0)/30,period[1/s]) 
A, Geometry 1 a , ™ 
#8 Materials 
Wy Heat Transfer (ht) 
a Heat Transfer in Solids 1 Entry method: 
“| Thermal Insulation 1 Start: 
=") Initial Values 1 ————— 
. > Heat Source stops pects) 
Ga Mesh 1 Number of values: 
ee Study 1 iE 
i, > Step 1: Time Dependent eee ey ee ee — q 














scl Configurations —— sail Gaiman 


12. | Under Time-Dependent Solver, make Relative tolerance to be 1.0e-5. Left-click 





on bead to compute. 





“Wl Medel Builder = = 1) [je [i Model Library) & Material Browser | =" | 
©@ 2-D.mph (root) —  — | Gl 
& Global Definitions ; | 
FP. Parameters Kk. Time-Dependent Solver 
Model 1 fred 7) 
= Definitions = General 
i Negara Tires: range(0,(period[1/s]-0)/30,pericd[1l 5 
wi Materials | 
SH Heat Transfer (ht) Relative tolerance: 106-5 
Heat Transfer in Solids 1 
|) Thermal Insulation 1 Absolute Tolerance 
=) Initial Values 1 » Tame St - 
™ ‘tepping 
) Heat Source 1 
Ge Mesh 1 - Result: While Sohing 
eat Stucky 1 
- Output 
i = Shep li Time Dependent {ou 
Pree coher Configurations . Adherent 
Fe] Solver 1 


24+ Compile Equations: Tin|| * boo 
iW Dependent Variables 1 
Ie. Time-Dependent Sole 
=] Direct 
hy Advanced 
ae Puolhy Coupled 1 
Direct 1 
foe Results 


14. Finally, under Results, select Surface under 2D Plot Group and select time to be 1 
only, we should be able to see the result like Figure 25. 
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® 2-D.mph - COMSOL Multiphysics 





File Edit Options Help 


0286 08 4- &. 


*oesgcue 








CA Report 














14. Further data analysis can be done under Report, such as generate a movie and export 
data and further compare results with MATLAB. 























0.8 


0.6 


0.4 


0.2 


(TT Model Builder \_ 7 = G1) {Settings > ll] Model Library| % Material Browser)  Gi|/¢h Graphics \__ 
(9 2-D.mph (root) Jf | a 
= Global Definitions aa 
P; Parameters (4 Surface 
() Model 1 (mod1) 
= Definitions — 
® Materials 
Fiettontene | we 
) Heat Transfer in Solids 1 
& ) Thermal Insulation 1 
v Expression el 
Initial Values 1 - sitai 
- \ > Heat Source 1 Expression: 
&3 Mesh1 T 
€3 Study1 
(0. > Step 1: Time Dependent Unit: 
fi] Solver1 mice 
sit Compile Equations: Tin eee 7 
ux Dependent Variables 1 Temperature | 
(iQ, Time-Dependent Solve 
ES Direct > Range 
‘y Advanced Ginae sia 
a= Fully Coupled 1 iii 
'S Direct 1 Coloring: Color table ’ 
{i Results 
#8 Data Sets —— 
s+ Views (V] Color legend 
8.85 
en Derived Values Fl Wireframe 
Ee Tables 
By 2D Plot Group 1 > Quality 
(©) Surface 
> Inherit Style 


-0.2 


QQ A6|m@° 5) 








0 
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0.2 


Surface: Temperature (K) 


0.4 0.6 


0.8 


1 


12 


ie 
a 343.4 


300 


200 


150 


100 


0 


v -4.975x10-? 
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APPENDIX G. COMSOL CODE FOR 3-D SIMULATION 


I, shames COMSOL 4.0a with 3D and hit = 





Select Space Dimension o> 


@ 3D 
. 2D axisymmetric 
2D 
1D axisymmetric 
810 
‘OD 
2. In Heat Transfer, select Heat Transfer in Solids (ht) and right-click Add selected then 
hit =. 
Add Physics 


> =, Recenthy Used 

> & AC/DC 

pe 2] Acoustics 

Ses Chemical Species Transport 
rs LI Electrochemistry 

t #2 Fluid Flow 

a (\\ Heat Transfer 





Heat Transfer in Solids (ht) 

BS Heat Transfer in Fluids ( oF Add Selected 
24 Heat Transfer in Porows lwreore - 
« Heat Transfer with Sire des Seataks: Radiation (ht) 
ail Bioheat Transfer (ht) 

—e- Surface-to-Surface Radiation (rad) 

a Joule Heating (jh) 
> Go Plasma 
& Au Mathematics 




















+ 
Selected physics 





=a Heat Transfer (ht) 
3. In preset Studies, select Time Dependent and hit Finish *. 


ie Model Library Fi 





Select Study Type do cb | | 
Studies 
@ 8) Preset Studies 
[= Stationary 
(. Time Dependent 
b> &) Custom Studies 





Selected physics 


| = Heat Transfer (ht) | 
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4. Under Model Builder, right-click on Global Definitions and left-click to select 
McLean ae Parameters as following: 
sel Bolder =F) rt Settings | Se El = 









r = 3-D meth. yee 


*», Parameters 
4 —& Global Definitions 





P| Parameters 

» WE Model 1 (mod) 

> fm Study 1 

> fi Results d 0.02[m] 002 m 
period L[s] ls 
a 0.25[rm] 025 m radius of rotatin...y in x direction 
b 0.25[m] 0.25 m radius of rotatin...y ony direction 
i) 0.5{m] 0.5m center of rotation in x direction 
yl 0.5[m] 0.5m center of rotation in y direction 
i) 2Oe5(W/m*2] 500000,0 Wy'mn* 


5. Right-click on Definitions and left-click to select Variables. Input Variables as 
following: 


‘Ty Model Builder ~~ 5 || (ie Settings a = 
(9 3-D_mesh.mph (root) 
= Olobal Definitions 


V) Model 1 (mod) Geometric Scope 
= Definitions 


sabes | Geometric entity level: Eniemodel 


ib Boundary Syste 





a= Variables 





ML View1 + Variables 
A, Geometry 1 
a Materials Name | Expression Unit | Des..ien 
® Heat Transfer (ht) XC wl+a"cos(2"pi"t/ period) m 
> Mesh 1 yc yO+b"sin(2"pi"t/ period) nm 
ee Study] q 1) (2"pi"d 2] fv 2])"expl-(tx-ne)} 2+ (y-ye) 22d“ 2) Wim! 
ie) Results 


6. oe aus on Geometry and left-click “Work Plane” and input | in z coordinate. 
| il me | || a2 Settings ~~ 


»“ Geometr 





a ©@ Untitled. mph (root) 
= Global Definitions 
@ &8 Model 1 modi) 
> = Definitions 


- Geometry S« 


ail | Geome ry l | — 
bias os ay) | Build Al i 
S@ Materials 


+e Import 
D ¥ Heat Tra 


3 Meshil *—) Block 
> G2 Study 1 +p Cone 
> Mil) Results * Cylinder 
* > Sphere 


More Prmitres La 


+S Work Plane 
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7. Right-click on Geometry under Wok Plane, select Circle. We need to build up two 
circles. One has radius 0.32m and the other one has radius 0.18m, both are centered at 
x=0.5m and y=0.5m. Right-click on Geometry and add Extrude, select Reverse direction 
and input distance to be 0.05m. Right-click on Geometry and add Block with Width, 


Depth and Height all equal 1. 
r FA Geometry 1 
a 46) > Work Plane 1 (wp) 
YA, Geometry, 


> Na View2 | [i] Build An 


> Form Union 4 a 








a2 Materials port 
3 Heat Transfer (ht) * CO) Circle 
<3 Mesh1 +o Ellipse 


. ee eed TV 
First circle: 
@ Untitled. mph jroot) 
= Global Definitions 
4 &%8 Modell (rod) 
= Definitions 
a ¥\) Geometry 1 
ai = > Work Plane l fwpt) 
al rs Geometry 
e] Circle d fet) 
kL} View 2 
wa > Form Union (fin) 
She: Miaterials 
= Heat Transfer (ht) 
SS Meshl 
am Stuchy 1 
Results 


Second circle: 
2 2 Untitled.rmnph (root) 
E Oolobal Definitions 
@ 8 Modell (rod) 
SF Definitions 
ai sy Geometry 1 
a 2 > Work Plane 1 fwpt) 
a Fs, Geometry 
Co) Circle 1 fet) 
Circle 2 fo) 
> bl View 2 
wa = Form Union (firy) 
Mee Mabterrnals 
= Heat Transfer fhe} 
fad Mesh 1 
> fm) Sturdy 1 


Extrude: 






























































C4 Circle 

~ Object Type 
~ Size 

Radims: O.32 


~~ Position 


Base: | Center 


re 0,5 
Wi 05 
¢ | Circle 


~~ Object Type 
Type | Solid 
~ Sree 

Radius 0.18 


ww Position 


Rage: | Center 
= 0.5 
yO O5. 


9] 


Model L fro l) =" 


= Definitions 
YS Geometry 1 Work plame: | wp 








warplanes wr?) || aparece 
a, = Form Union ffir) | 
Ee Materials 
BH Heat Transter phe) [4] Keep input objects 
foe Mes 1 (| Keep cross-sectional faces 
Stuety 1 i | 
fea Results ~ Distances from Work Plane 





= 


Reverse direction 


Add a Block: 


“PD Model Builder ~~ _ ~ = E484 Settings ~~ __ 





= sppendim g.mph frect) 
Global Vetinitions a3 Block 
Kaeo! 1 frrererd ID SS 
=> Definitions 
PM, Gaccorrebiry 1 TrPpe [Solid , 
a Work Plane il (wept) | 
GE Extrude 1 fexta) w Ste: creed SP ncaprer 
ij Block 1 (fbikt) Wrielth: 1 
a, = Form Union ffir) 
6 Moaterial= cae | 2 
aay Stuchy 1 - Position 
Results wz , 
Aue hype | Cartesian 
bt Li 
We o 
zz iL 


FF Fastiiaticom Sigh 
e Layers 


The geometry should be built like this: 
POG (Aan etek ek b|e seo 2 | ee) ee 
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7. Right-click on Materials and left-click on add Materials. Add Steel AISI 4340. 
Right-click on the geometry in Graphics and left-click on it to select. Make sure that all 
the portions are selected. “ 1” , “2” and “3” should be under the selection. 


OTT Model Builder * 

2 0) 3-Domech.mph fre) 
¢ ££ Global Definitiers 
ao) Baodel ll feel 

. © Definitions 
. A, Geometry] 
a 2 Moatenal: 





| & Maternal 


Geometric Some 
Geometric entity level Domain 


Selector: 


» BE Sheed AMS 434d) 
Heat Transfer fide! 
+ a Kaew J 


» B® Derveed Values 
» ES Taller: = 
» Ti 30 Pint Greap i 
+ fm 3D Piet Greop s 
» 9 Probe 1D Plot Grengp 5 


All corresens 





F Raaterial Properties 


37 Matertal Contents 


a El Pace ceit aad _ aa : — 
i! Plaverl Property | Marne | Value | Wrirt | 


ia) Gnrrestond 


























Thermal comduictety k 


475[.7K)) J fea) 
rh PAD) kg. "3 
44,5,,.°K}] eee 


8. Right-click on Heat Transfer and left-click to add Boundary Heat Source. For 


6G, 99 


Boundary heat source Qp, enter “q’’, 
Right click on selection 9 on Graphics and left click on it to select. 





(8 3-D_mesh.mph (root) 
= Global Definitions 
4 4) Modell (modi) 
= Definitions 
» YA Geometry 1 
Materials 
4 8 Heat Transfer (ht) 
> Heat Transfer in Solids 1 
% | Thermal Insulation 1 
S” Initial Values 1 
_) Boundary Heat Source 1 
G3 Mesh1 


» BS Derived Values 
» HE Tables 


which is the heat source defined from the Variable. 


_) Boundary Heat 


Boundaries 
Selection: | Manual 
q 


* Boundary Heat Source 


Boundary heat source: 


On gq 
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1] eh Grephics 


| Ge #-)QQAeBl>-he =| 2 oer 





) Boundary Heat Source 
Boundaries 


Selection: | Marat 





* 
z ty te 
t - 
n x 

Rounclans heat srr 
Gn q wie 


in the beginning. 
a OW Heat Transfer (he) 
“> Heat Transfer in Solids 1 
| Thermal Insulation 1 & x 
© Initial Values L 
) Boundary Heat Source 1 
t Gat Mesh 1 
t ee Shody 1 
a (ii Results + Initial Values 
» #! Data Sets 
» 2% Derived Values Temperature: 
> SE Tables T OK] 
> (i 3D Plot Group 1 | | 
> (i 3D Plot Group 2 ee 
> PX Probe 1D Plot Group 3 10 


Wii? 
10. Right-click on Mesh and select Free Tetrahedral, repeat this process three times. 


a 
oe 
of 


r= 
[ 








More Operations 


First Free Tetrahedral mesh: selection 2 1s picked and makes the predefined Element Size 
to be Extremely Fine: 
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tty Madd Bailes 7 © Game 
a 81-1) mest. rea 
# Geta Beaice: 
al i Ficoel L freon’ li 
» = Ushrihore 
A, Goudy) 
, hee 
pW Hest tovete: rc 
i i Meal 


BBE AO) bape. 
POé@e- TR AeRM dei le |e fier 





9 SE eraeel alee 
 Teble 
» IDA tral 
@ Wana 
= Preke ID Pit Gerd 
Fy Pl bapst 
= Peal 
a Auruaiers 


@ Predetred |Ectsurnaby firm 

















Second Free Tetrahedral mesh: selection 3 is picked and makes the predefined Element 
Size to be Extremely Coarse: 


i five ela. — 


a bed Ld ree'l pees a ’ 
= Dehnitione a 
Fs, Genesetiy | Geortne at: Heat | Cerny sm] 
® Mae 


Eicrecal Sur 


[El Teme: ©) Caran 
b ikonend Soc farecton 
i) Free 10 Four fevagn 3 
a OF) Peoot. 
(2) Peper 
EE ftatinn | 


Third Free Tetrahedral mesh: selection 1 is picked and makes the predefined Element 
Size to be Extremely Coarse. 


Sineta Niualer = 


Bite 
EGG *- O00 448 \+-s ei) a ie, Bee 





es Bi pashungh fred a = 
© Cabal Petates: a-are 
a (bE Badd] jeer! 


; 4 Pret Teed i! 
a Gen! 

ijn 

& Tebcheael 





Si Prete De Sed 
a 2) Raat 
be Phar 
Ee Serkierl 








a to build the mesh. 
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Then select build all button 


11. | Under Study in Step, select Range button ' , under Entry method, select 
Number of values, start from 0 and stop at 1. Put Number of values to be 51. 


08 3-D_mesh.ngh (root) _ 
= Global Definitions i. Time Dependent 


v Model 1 (rod) ~ Study Settings 
f= Study 
\(\, > Step 1: Time Dependent Times: —_range(0)1/30,1) 5 thi 
Tr, Sohver Configurations 
fai Results 
#2 Data Sets 
EE Derived Values 
El Tables 
fl 3D Plot Group 1 
(i 3D Plot Group 2 
fi, Probe 1D Plot Group 3 
EE Report 
= Playerd 
@8) Animation 1 





12. Under Time-Dependent Solver, make Relative tolerance to be 1.0e-5. Left-click 


on to compute. This 3D problem may take couple hours to solve, in order to see a 
quick solution, we can make all meshes to be extremely coarse. 



































"| Medel Builder =i | <Q 
3-D_mesh.mph {root} = 
” = Global moh ks Time Dependent Solver 

&) Model 1 fod) —= = 

i Stucky 1 Gener i = 
[. = Shep 2: Time Dependent Tirrnes= earige(o,1/50,1) bs [tat | 
he. Sobver Configurations 
= as " Relates tolerance le-5 


Sc Compile Equations: Time Dependent 


F Absolute Tolerance 
ae Dependent Variables 1 


Se liene-Dependent Sober] b Tome Steppang 
Time-Dependent Solver 3 Tame Steppe 
a) Direct 
% 5 seat b Results While Solving 
atm Fully Coupled 1 » Output 
(S] berate l 
ig Results F Adve eel 
ED Data Sets = 
HH = 


BA! Oerved Valucs 


2. Finally, under Results, select Surface under 3D Plot Group and select time to be 1 
only, we should be able to see the result like Figure 28. 


Meddbeko + Satie 1 ~ Ole 


8 340 raturegh (eed 
ss Eiloied issami AD Plot Group 
Gi Mid | fered] = Dee 
me Tink | = 
@ toe det: take Ce 
= Uris Sate ; ; 
ee Deed Nika Tite: I 7 
i Tabby 
Ce DE Pet ep > Ae Sleep 
2 ten — 
Pet Gee? ‘ear fedora rrrrrr— 
Pe crepes fate 7 _ rt ax =~ |. pwedla 
Lo) Figs ; 
i Seercril F slse | ick z 


Fh here [lene bik 4 * 

















Phe be Wa hee: Sere 





14. Further data analysis can be done under Report, such as generate a movie and export 
data. We can add a Domain Probe Point under Definition to analyze the temperature rise 


EN [aia . PR TD dee FO | Ae) Lee @ ofa 

















2 G8 Ll vetorgi eat 
a F Sl doee 
8 Beery 
ot HD irked Paar Li 
aE Defesiccs 
= Senekemt 
ie Rades Spas hp 
b oft ime 
eS, beech 
) Bod 
p Hear Teece pa 
bp ie Baek. L 
ip Shey 
Gia Aer lly 





RES : ———— 
Click on Probe 1D Plot Group under Results, we can plot temperature as a function of 


ee 


. oe? | QL = Spt congener aq nei iille=t 
Shei |S. fable Pint | Fama pat Mi yo rote cage 





Tareas hi 
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